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NOMENCLATURE 


A - Combustion chamber cross-sectional area 

a - Downstream nozzle area 

Ci,C 2 ,Cu - k-e model constants 

0 - Combustion chamber diameter 

d - Inlet nozzle diameter 

k - Kinetic energy of turbulence 

L - Contraction nozzle downstream distance 

z - Length scale 

R - Inlet nozzle radius = d/2 

s - Swirl number 

v=u,v,w - time-mean velocity (in x-,r-,9 directions) 
x,r,9 - Axial, radial azimuthal cylinderical polar coordinates 

X - Downstream distance from domain inlet 

e - Turbulence energy dissipation rate 

* - Swirl vane angle with respect to facility axis (except in 

Chapter III) 

u - Turbulent viscosity 

p - Time-mean density 

¥ - Stream function 

to - Vorticity 

y - Kinematic viscosity 


Subscripts 

n - Relating to Swirler hub 

m - Station maximum value 

o - Value at inlet to flowfield 

x - Relating to axial direction 

0 - Relating to swirl direction 


Superscripts 

( ), - Time-mean average 

( ) - Fluctuating quantity 


iv 


CHAPTER I 


INTRODUCTION 

1.1 Turbulent Swirling Flow and Its Observed 
Physical Features 

Turbulent swirling flows are of considerable practical importance. 

In heat generators, the designers of combustion chambers and industrial 
furnaces generate swirl to assist in the stabilization of the flame and 
to promote rapid mixing. In heat transfer devices, the swirling flow is 
utilized to increase the rate of heat transfer. And in nature, swirl 
creates a drag force by causing trailing vortices which occur when the 
separated flow over the aircraft wings rolls up into vortex cores at the 
wing tips. The swirl can be imparted to the flow by different means. 

The physical geometry and aerodynamic purpose are the main factors which 
lead to the method of imposing swirl on flowfields, adjustable vane 
swlrlers, rotating the inner or the outer cylinders of annular flows, 
tangential inlet jets, and other techniques. The degree of swirl 
usually is. characterized by the swirl number S, which is a nondimen- 
sional number representing axial flux of swirl momentum divided by axial 
flux of axial momentum times equivalent nozzle radius. It is described 
at length in References 1 through 3 and in Section 8. 2. 2. 2 of Chapter 8. 

Swirling flows have been classified by King et al . (4) according to 
the wall boundary conditions as follows: 
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1. Unconfined swirling flows where the wall effects are 
negl igible. 

2. Confined swirling flows in short, large-diameter chambers where 
side wall effects strongly interact with the swirl to produce signifi- 
cant secondary flows. 

3. Confined swirling flows in tubes where circumferential wall 
effects interact strongly with the swirl flow. 

The time-mean and the turbulence characteristics of swirling flow 
are generally affected by the wall boundary. A swirling flow in a 
stationary duct or annulus possesses the characteristics of a vortex 
motion with axial motion along the duct or annulus. The swirl motion 
decays, and the flow characteristics approach those of a pure axial 
flow. The axial and swirling velocities of swirling flow in an annulus 
with the inner cylinder rotated, develop with distance along the annulus 
and become fully developed at certain lengths (5). 

Vortex breakdown has been observed in the swirling flows through 
nozzles and diffusers, and in combustion chambers. The occurrence and 
place of vortex breakdown depend on a balance between the magnitude of 
the swirl, the external pressure gradient, and the degree of divergence 
of the flow. The greater the adverse pressure gradient, or the degree 
of divergence, the lesser is the swirl that is needed for the vortex 
core to breakdown, and the higher the swirl, the more the breakdown 
moves up into position from downstream (6). A vortex core forms after 
the occurence of the vortex breakdown phenomena and establishment of a 
reverse flow zone. It processes about the axis of the flow on the 
boundary between the zero velocity and the zero streamline. A variety 
of forms of vortex breakdown and their characteristics has been recently 
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identified by Leibovich (7). Most recently, a book edited by Swift 
et al. (8) has brought together recent developments in studies of vortex 
flows. 

Experimental and numerical studies (2, 9-11) show that swirl has 
effects on turbulent flows; jet growth, entrainment, and decay; and 
flame size, shape, stability, and combustion intensity are affected by 
the degree of swirl imparted to the flow. It is found that the swirling 
jets spread more quickly and the mean velocity decays more rapidly than 
in a nonswirling jet. In duct flows, at high swirl strength, back flow 
can occur in the central region of the flow. The central back flow 
plays an important role in flame stabilization by providing a hot flow 
of recirculated combustion products and a reduced velocity region where 
flame speed and flow velocity can be matched. 

1.2 Practical Need and Method of Investigation 
.of Swirling Flows 

Because of the inherent importance of swirling flows in many 
engineering applications, their investigation has been the subject of 
continued interest of many fluid mechanic researchers. The field of 
turbomachinery is among the various implementations of the swirling 
flow. The increased interest in high performance of this equipment 
produced the need to minimize the losses through the engine components. 
Chief among these components are the combustors. Combustor losses are 
mainly attributed to pressure drop, incomplete combustion and flame 
stabilization problems. The most fruitful avenue of approach to 
eliminate these losses is through understanding the phenomena which 
generate them. This can be accomplished through accurate knowledge of 
the characteristics of the combustor design parameters. 
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The understanding of the detailed nature of the behavior of the 
swirling flows can be aided not only by experiments, but also by their 
predictions. Early attempts to predict swirling flows with simple 
boundary and geometric conditions, and with constant fluid properties 
followed the integral approach, that is, velocity profiles were assumed 
for the flow and the integral equations of motion were solved (12, 13). 
Recently, a more general prediction method has been used. This solves 
the motion governing differential equations rather than the integral 
equations via mathematical models incorporating numerical finite differ- 
ence computer codes. Improvement and use of these techniques will 
significantly increase understanding and reduce the time and cost of 
development. 

The present study forms a theoretical part of a project on the 
investigation of flowfields found in typical combustor geometries. The 
partial differential equations of conservation of mass, momentum (in x, 
r, and 0 directions), turbulence kinetic energy and its dissipation are 
solved using an advanced version (14), called STARPIC, of the Imperial 
College TEACH-T computer program (15). This code solves the correspond- 
ing finite difference equations of the motion governing equations using 
a semi-implicit line-by-line method for values at points of a variable 
size rectangular grid, with variable under-relaxation. The turbulence 
characteristics are described by a modified version of the k-e turbu- 
lence model. This code will be described in more detail in Chapter 3. 

1.3 The Present Study Objectives 

The present study is directed towards an understanding of the aero- 
dynamics of the flow occuring in an isothermal can-type combustor 
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flowfield. To attain this goal of more accurate prediction of the 
characteristics of the flowfield under consideration, the following 
objectives are considered. 

1. Evaluating the sensitivity of the predicted results to dif- 
ferent inlet boundary condition types. 

2. Investigate the effects of the degree of contraction of the 
flow at the exit of the domain. 

3. Develop an advanced version of the k-e turbulence model to more 
accurately predict confined turbulent swirling recirculating flows. 

4. Perform experimental tests as needed to verify the precision of 
the theoretical developments. 

The first of these is intended to demonstrate the need to use 
realistic inlet boundary conditions in order to predict more correctly 
the subsequent flowfield. The second objective increases the range of 
application of the code by the inclusion of downstream contraction 
nozzle effects. The development of a suitable advanced turbulence model 
from detailed turbulence measurements is the subject of objective 3. 

The resulting powerful modal is to be expected to simulate the flowfield 
more realistically, and predictions are made to demonstrate this. 

From an economical point of view, it is necessary to use compu- 
tational fluid dynamics procedures rather than experimental approaches 
for a better definition and understanding of the flowfields in the 
combustor. However, the most viable approach would be for a parallel, 
interactive, theoretical and experimental program to be used. Thus, the 
last objective of this study is to perform experimental tests as needed 
to verify the precision of the theoretical developments. 


1.4 Outline of the Thesis 


The first chapter of this nine-chapter thesis is the introduction. 
Turbulent swirling flows and their observed physical features are 
briefly described. This chapter highlights some of the practical needs 
and methods of investigation of the swirling flows. Finally, the objec- 
tives of the present study was stated an# justified. 

Chapter II presents the method of investigation of swirling flows. 
Both numerical and experimental studies of other workers are explored. 
The experimental study findings and the accuracy of the predictions of 
those workers ire discussed. 

The flow governing equations, the grid system and the computational 
techniques employes? in computer code (STARPIC) of the present study are 
presented in Chapter III. The boundary conditions imposed and the code 
deck are also discussed in this chapter. 

Chapter IV provides a brief description of the structure of turbu- 
lent flows and the mathematical models which are expected to represent 
the flow behavior. This chapter is concerned with the survey of the 
contributions of many workers to the understanding of the characteris- 
tics of turbulent flows, and with the turbulent flow structure and its 
modeling approaches. Various turbulence models are discussed: Prandtl 

mixing length, energy, one- and two-equation models, algebraic and 
Reynold stress closures. Among the two-equation models, is the k-e 
turbulence model. This model is the subject of the following chapter. 

The structure of the k-e turbulence model and method of evaluating 
its constants are reported in Chapter V. This chapter exhibits the 
previously proposed extensions of the k-e model and the limitation of 
each extension. Chapter V provides the background for the extension of 
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the k-e model to applications in turbulent confined swirling flows, as 
shown in Chapter VII. 

Chapter VI, on the other hand, analyzes recent turbulence measure- 
ments, provided by another worker, to investigate the flow boundary and 
swirl strength effects on the characteristics of the turbulence param- 
eters. Also this chapter provides useful information that helps in 
extending the k-e model to the cited flow type. 

A validation study of the previously proposed extension of the k-e 
turbulence model to swirling flow applications is conducted in Chapter 
VIU As a consequence of this validation, two extensions to this model 
are proposed: the first is accomplished by optimizing the model con- 

stants according to the available data; the second provides an empirical 
formulation for the k-e model constant Cy The predictive capability of 
both extensions are tested by comparing the results with measurements 
obtained by other workers using a 5-hole pitot probe and a hot wire 
anemometer. 

Finally, Chapters VIII and IX are the closing chapters of this 
effort. The former presents general predictions; the effects of imposed 
boundaries and swirling flow strength on the flowfield characteristics 
are all investigated at length. The boundary types include the inlet 
flow profile specifications, and the expansion and the contraction 
ratios at the inlet and at the exit of the flow domain, respectively. 

The latter. Chapter IX, summarizes the main conclusions of the present 
study. 

There are also two appendices which provide the tables and figures 
which are presented in this study. 


7 


CHAPTER II 


PREVIOUS STUDIES IN TURBULENT SWIRLING FLOW 
2.1 Introduction 

Numerous theoretical and experimental swirl flow combustion studies 
currently exist. The former have been obtained through the use of math- 
ematical modeling of the combustor flow domain; usually finite differ- 
ence techniques are employed in such modeling processes. Several text 
books (16, 17, 18), recent papers (11, 20-22) and more recently an 
international conference (19) exemplify studies and give soluable 
insight into our understanding of the swirl phenomenon. On the other 
hand, the latter have been conducted using different measurement tech- 
niques. Chief among these techniques are pitot probes, hot-wires and 
laser anemometers. 

2.2 Numerical Analysis of Swirling Flows 
2,2,1 Numerical Techniques 

The flow of an incompressible viscous fluid in two space dimensions 
is most often treated by the stream function-vorticity tp-m formulation. 
However, the grid nonuniformal ities and high swirl produce problems in 
this technique (22). Also, the use of the stream function and vorticity 
as the dependent variables causes complications and difficulties in 
implementing the boundary conditions for closed systems. The advantage 
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of the stream function-vortici ty technique in reducing the number of 
partial differential equations may be outweighted by the above 
difficulties. 

To overcome the problems mentioned above in using 'J'-oj formulation, 
more recent emphasis is being placed on direct solution for primitive 
variables p, u, v, and w instead of t-ui. This is based on the numerical 
solution of a finite difference representation of the Navier Stokes 
equations in primitive variable form. 

Examples of the previous studies which have employed the afore- 
mentioned techniques are summarized and displayed in Table I. The 
essential differences between the various computer codes which were 
employed in these studies, are the method of deriving the finite differ- 
ence equation for the set of governing equations, which were selected to 
simulate the physical process, and the location of variables in the grid 
system, 

2,2,2 Accuracy of the Predictions 

The capability of current computational fluid dynamic codes has 
been evaluated by many workers. Most recently, Sturgess (40) evaluated 
the physical modeling embodied in such codes. Such evaluation was 
achieved by a detailed comparison against related experimental data. It 
has been demonstrated that qualitatively correct predictions have been 
obtained for a variety of complex flows. Absolute accuracy has not yet 
been demonstrated. The discrepancy may be attributed to the physical 
model and numerical methods employed. Following are samples of the 
recent validation work found in the literature. 
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Hendricks and Brighton (41) investigated mixing a turbulent jet 
with a coaxial slower-moving secondary stream in a constant diameter 
tube. The analysis involved a numerical solution of the governing flow 
equations which were simplified by the Prandtl boundary layer assump- 
tions. The effective viscosity was calculated from a modified version 
of the two-equation turbulence model. They found that the predicted 
centerline velocity lags the experimental data by about one pipe radius 
while the maximum pressure point is arrived at about three radii too 
soon. However, the value of the predicted maximum wall pressure is in 
excellent agreement with the data. Also, it was found that the predic- 
tions of swirling flow were more accurate than of nonswtrling ones. 
Generally, good agreement between the calculated results and the corre- 
sponding experimental data was obtained. 

Predictions of isothermal and reacting flows in two dimensional 
combustor geometries were reported by Jones and McGin'rk (42). Turbulent 
transport was approximated via a variable density form of the k-e 
turbulence model. Their predictions for a two dimensional gaseous- 
fueled combustor indicated that the physical models used were capable of 
calculating the flow oattern and heat release fairly well. 

An evaluation of currently available turbulence and combustion 
models as applied in gas turbine-like flows was performed by Syed and 
Sturgess (43). They used the k-e model and the eddy breakup model in 
their predictions. They found that the turbulence model predicted mean 
flow quantities to an accuracy of about 85 percent in regions of para- 
bolic flow and qualitatively wel 1 -predi cted in elleptic flow regions. 

The size and strength of the recirculation zone were underpredicted by 
about 20 percent. 
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Morse (44) conducted theoretical and experimental work in order to 

'i 

i 

evaluate the performance of turbulence models as applied to turbulent 
swirling flow. Performances of mean flow closures (Reynolds stresses 
components are assumed to be locally determined) and Reynold stresses 
closure (Reynolds stresses are obtained from transport equations) were 
discussed in detail. He reported that the difference between the rate 
of spread of a round jet and a plane jet issuing into stagnant surround- 
ings cannot be predicted using mean flow closures. In addition, these 
models fail to reflect the experimental observed increase in spreading 
rate with swirl number. Also, he concluded that the Reynolds stresses 
model failed to bring good agreement for corresponding olane and 
axi'j-metric shear flows. I 

A numerical simulation of the isothermal flow in a diffusion-flame f 

research combustor was carried out by Sturgess and Syed (45). An 
axisymetric finite difference solution of the time-averaged, steady- 
state, elliptic form of the Reynolds equations was used with closure 
provided by a two-equation turbulence model. This simulation gave 
acceptable agreement with experimental data in many respects, except the 
prediction of the recovery rate of a central near-wake region. The 
authors provided detailed explanation of the source of the discrepancy. 

Generally, the findings of these studies show that turbulence 
modeling and numerical errors are not great enough to allow false con- 
clusions to be drawn about the main physical characteristics of the 
flow. Also, it should be mentioned here that turbulent swirling flows 
are complex, and the level of accuracy and reliability of the predicted 
results strongly depend on the complexity of the flow. 

■ . i. 

. . i 
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2.3 Experimental Investigation of Swirling Flows 


2.3.1 Measurement Techniques 

Experimental observations of swirling flows have been made for a 
number of years. Their preliminary approach was by means of a pitot- 
static probe which consists of an open-ended tube placed to face the 
direction of the flow, together with another open-ended piece of tube 
placed at right angles to the direction of flow. Later, more elaborate 
equipment was used to measure the magnitude and direction of fluid 
velocity. This instrument is the 5-hole pitot probe. As more detailed 
flow characteristics were needed, the hot-wire anemometer was used. 

Very recently, Perry (46) described this tool in detail, in terms of a 
variety of constructions of hot wire applications and theories. 

Finally, the advent of a laser doppler velocimeter (LDV) has had a 
significant impact on data acquisition in complex flow under the combus- 
tion environment. Table II displays a collection of the previous 
studies which employed the aforementioned measurement techniques. 

2.3.2 Experimental Study Findings 

The relevance of the swirling flows to the aircraft industry has 
encouraged many researchers to study its characteristics. Mostly, the 
findings of these studies have been used for the evaluation of the 
capability of current computational fluid dynamics codes for use as 
a design and analysis tool in gas turbine combustor development. 

The following selected studies are provided to exemplify such 
investigations. 
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Murakami et al . (47) performed an experimental investigation of 
swirling flows in cylindrical ducts. Two different types of swirl 
flows, namely a forced and a free vortex type, were introduced in the 
test section. The former was introduced by an axial flow type impeller 
located in the swirl chamber, and the latter was developed by a radial 
guide vane swirler. It was found that the cited swirling flow shows 
three kinds of flow patterns in its course of decay: forced-free vortex 

type, transitional type, and forced vortex type. Also, it was found 
that the swirl strength decreases exponentially along the pipe axis and 
its decrement varies with the change of the flow pattern. 

Swirling flows in stationary channels were investigated by 
Sukhovich (34). A stream of air was supplied to the working zone by a 
nozzle of area 0.7 times the channel area. The secondary, swirling flow 
was imparted to the flow through an annular gap, formed by the outer 
surface of the nozzle and the inner surface of the working tube. His 
analysis of the experimental data showed that the axial and tangential 
stresses are anistropic and the turbulent viscosities depend on the 
swirl strength and coordinants. Also, he found that at high gradient of 
annular velocity and Richardson number, the turbulent viscosity for the 
tangential component of momentum is lower than the corresponding one for 
axial component of momentum. 

Measurements in a model combustor composed of two confined coaxial 
swirling jets under noncombustion conditions were reported by Vu and 
Gauldin (38). Swirl for the inner jet was generated by a small axial 
flow swirl generator with fixed vanes, and that for the outer jet was 
imposed on the flow by adjustable vanes. Detailed time-mean and 
fluctuating flow measurements were obtained for coswirl and counterswirl 
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conditions with a directional pitot probe and hot-wire anemometer. It 
was noticed that the outer swirl has a strong effect on the formation of 
a recirculation zone and on mixing characteristics in the inner jet 
shear layer. Also, they found that high levels of turbulent fluctua- 
tions and large dissipation rates characterize the central flow region 
for both coswirl and counterswirl conditions, while the outer flow 
regions exhibit relatively low turbulence levels. 

Confined turbulent swirling flow data obtained from a single hot- 
wire using the six-orientation techniques was reported by Jackson and 
Lilley (10). Swirl was Imparted to the incoming flow by means of a 
variable angle flat vane swirl er (3). The effects of swirl strength and 
degree of contraction of the flow at the exit of the flow domain were 
carefully examined. The effect of swirl on the time mean velocity field 
was found to shorten the corner recirculation zone length and to 
generate the existance of a central recirculation zone, which is 
followed by a processing vortex core region. The central recirculation 
zone was reduced and the turbulence levels and shear stresses were found 
to increase by introducing a nozzle of area reduction ratio of 4 located 
2 chamber diameters downstream of the flowfield. 

2.4 Summary of the Chapter 

The material which is presented in this chapter exhibits the most 
relevant work to the subject of confined turbulent swirling flows. It 
provides valuable insight into the directions to be taken in the inves- 
tigation of such flows, and affords an up-to-date data base to help in 
validating and improving of the physical models and numerical 
techniques. 
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CHAPTER III 


THE STARPIC COMPUTER CODE 
3.1 Introduction 

In recent years, there has been considerable progress in the com- 
bustor flowfield predictions using finite-difference numerical proce- 
dures. Several fluid flow computational techniques have been employed 
to improve the flowfield simulation and turbulence modeling. These 
techniques included method of approach, grid system, numerical stability 
and many others. 

The computer program STARPIC (acronym for Swirling Turbulent 
Axisymetric Recirculating flow in Practical Isothermal Combustor geome- 
tries) is among the previous codes which are displayed in Table I. Most 
recently, the code has been expanded to include reacting flows (29, 48). 
The isothermal version of this code has been selected to be employed in 
the theoretical aspects of the present investigations. This selection 
was made because STARPIC is the best code yet available for predicting 
isothermal flowfields in axisymetric combustor geometries. This code is 
an advanced version of the TEACH computer code (15). A complete 
description of the finally developed computer program, with a full 
explanation of the equations, source terms, revised cell volumes for 
axial and radial velocities, constants occurring and techniques for 
handling turbulent swirling flow near curved boundaries, is available 
(14). Therefore, it will be described only briefly In the following 
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sections, so that it will be available whenever is needed to refer to 
it. 

3.2 Flow Equations and Grid System 

The geometry under consideration and the corresponding experimental 
and solution domains are shown in Figures 1 and 2, respectively. It can 
be easily foreseen that recirculation zones will exist within the solu- 
tion domain and as a consequence, the governing equations are elliptic. 
The set of axisymetric governing equations may be compactly represented 
in the following appropriate form: 

It (o'* -r, |$) ♦ Ip [pvt - r + fi] - sj 

where 4> is any of the flow variables (u, v,...etc). is the corre- 
sponding effective exchange coefficient and is a source term. This 
equation may formally be integrated, over a typical cell, to yield: 

A* - * A"* - r + f ]dA = p I v s+dv 

where N, S, E, W and P refer to north, south, east, and west of the cell 
and the reference location in the cell, and A and V are the area of the 
cell boundary and the cell volume, respectively. If the variables are 
assumed to be linear inside the cell, the above equation can be 
rewritten as: 

a e - to* - Ify A w + f \ A M 

- [pv* - r t |i] s a s » [s** p * sj] x v 

where r^, Sp’ , and are tabulated for each 4) in Table III. The left- 
hand side of this equation represents the total transport across the 
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cell boundaries by convection and diffusion. This resultant difference 
equation has the following form (49): 

(A p - S p ) + p - E A n t n + S u 
n 

where: z denotes summation over neighboring nodes of a typical grid 

n 

node P; A n are influence coefficients; and Ap = z A n . The pressure is 

n 

obtained by combining the continuity and momentum equations in the 
manner explained by Caretto et al. (50). 

The flow domain is overlaid with a grid whose nodes are formed by 
the intersections of arbitrarily spaced coordinate lines, as shown in 
Figure 3. The u and v velocities, shown as arrows in the diagram, are 
located mid-way between the grid intersections. All other variables are 
located at the grid nodes. This staggered-grid system was originated by 
Harlow and Welch (51) and was earlier adopted by Patankar and Spalding 
(52) and many others. This grid arrangement has two special merits: 
first it places u and v velocities between the pressure that drives 
them, and it is easy to calculate the pressure gradients that effect 
them; and second, these velocities are directly available for the 
calculations of the convective fluxes across the boundaries of the 
control volumes surrounding the grid nodes (20). It should be noted 
that the axial velocity reference location is (1,0) even though it 
actually represents the velocity positioned at (1-1/2, J) 

3.3 Computational Techniques 

Solution of the system of finite difference transport equations 
with appropriate boundary conditions is achieved by means of an itera- 
tion procedure. Firstly an initial or intermediate value of the 
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pressure field is used to solve for an intermediate velocity field. 

Then, continuity is enforced, by solving the equations for pressure 
change and thereby determining the required adjustments to the pres- 
sures. Corrections to the pressure and velocity field are applied and 
the equations for the remaining variables are then solved in turn. The 
iteration is terminated when the correct solution satisfies the 
difference equations within a certain percentage of a global measure, 
such as the total mean flow rate or inlet momentum, depending on the 
variable in question (0.9 percent or less has been adopted for this 
study). 

The value of <J> at each grid node is related to its neighbors via 
i’p ~ Ajjj + A<- + [A^ 'J’e + ^ 

where N, S, E, and W are neighboring locations of P, the A's are 
coefficients combining the convection and diffusion effects of the flow 
and incorporates any existing source terms, which are calculated from 
the geometry, properties and dependent variables of the previous itera- 
tions. Convective and diffusive fluxes are determined from a hybrid 
difference scheme. Thus diffusion is neglected when the flow is highly 
convective and "upward difference" formulas are employed to express the 
differentials. On the other hand, when the diffusion is comparable with 
convection, the former is taken into account by employing a central 
difference numerical scheme. With the aid of a tri-diagonal algorithms 
and with the assumption that the quantity in brackets is known, all 4 > p 's 
may be calculated along any grid line. In this manner, one can traverse 
along all the lines in a vertical direction sequentially from left to 
right of the integration domain. 
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In order to avoid numerical instabilities, under-relaxation has 
been employed, that is 

* = f+j + 

where f is the under-relaxation factor 0 < f < 1 and f equals 1 when 
applied to the pressure field, <j> is the resultant value obtained from 
the present iteration cycle <j> j and the previous iteration cycle values 
♦j-l* The optimum value of this factor depends on the geometry and the 
boundary conditions under consideration. 

3.4 Boundary Conditions 

At the inlet plane, all variables are given definite fixed values. 
The axial and radial velocity profiles can be flat or an experimental 
distribution. The swirl velocity profile can be assumed to be that of 
solid body rotation or any specific distribution. The swirl number and 
the velocity profile types are the important features of the inlet 
velocity profiles. The value of k is assumed to be: 

k.j = 0.03 x (u 2 + v 2 + w 2 } 

where i is the grid location across the inlet, while e is calculated 
from: 

“t - c „ p X 2 /* 

After is estimated by mixing length theory as shown below: 

= k 3/2 /(0.005R) 
i 1 

■where R is the combustion chamber radius. 
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Near the wall boundaries, tangential velocities are connected with 
their zero wall values by way of the tangential shear stress wall 
functions, and the normal velocities are given zero normal ',radients. 

The e values are fixed using the length scale near the wall and the 
current value of k. The k and pressure are given zero normal gradients 
c^jrywhere except at the inlet. 

At the center line, the axial velocity and energy dissipation rate 
gradients, and the radial velocity are set equal to zero, while the 
swirl velocity is assumed to be that of the solid body rotation. The 
axial velocities at the outlet are deduced from their upstream values by 
adding a fixed amount which ensures overall mass conservation. The 
radial velocities are assumed to be zero while the other variables are 
assigned zero gradients. 

Since the flowfield boundaries are positioned midway between grid 
points, they coincide with normal velocities. This simplifies the 
insertion of specific boundary conditions by controlling the value of 
the coupling coefficients and the false source term. This will be clear 
if the difference equation is rewritten as: 

E a n " *n> = VP * S u 
n 

n = N, S, E, and W 

which implies that if any link, say aj is broken and Sp is set equal to 
zero, then S u represents the resulting flux of the neighbors. On the 
other hand, if the right hand side of the above equation is set equal to 



the specified boundary value. 
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3.5 Code Deck Descriptions 


STARPIC has thirteen routines in addition to the main code. A 
detailed description of these routines is available (14). Here a brief 
summary of the main features of each code are given. 

The main routine specifies the grid size and coordinates in addi- 
tion to the values of the field parameters at the inlet plane. The 
present study added a section which automatically generates a stepped 
grid for simulation of the exit nozzle. Once the grid size and coordi- 
nates and the domain boundaries are specified, subroutine INIT is called 
to set values of numerous geometric quantites concerned with the grid 
structure, and to initialize most variables to zero or other reference 
values. Also, the main code controls the flow to the other routines for 
each iteration as shown in Figure 4. After convergence is obtained, the 
final results are printed out and/or stored on magnetic tapes or disks. 

Subroutines CALCU, CALCV, CALCP, CALCW, CALCTE and CALCED are 
called in sequence for each iteration to update the axial and radial 
velocity fields, pressure, swirl velocity, kinetic energy and energy 
dissipation fields, respectively. Each of these routines calls PROMOD 
to modify the finite difference equation coefficients or other variables 
near walls or other boundaries where particular conditions apply. Once 
the call returns from PROMOD, routine LISOLV is called to update a 
particular variable over the entire flowfield by applying a TDMA tech- 
nique to all the lines in the r-di recti on sequentially from left to 
right of the integration domain. Finally, fluid properties are updated 
by calling PROP, and the totals of the absolute residuals are calcu- 
lated. This cycle is repeated until convergence is achieved. 
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3.6 Closure 


This chapter has presented the numerical techniques and solution 
procedure incorporated in the STARPIC computer program by v lich the 
governing difference equations of the flowfield are solved. The tech- 
niques employed are general and flexible and may be applied to calculate 
various flow situations governed by different boundary conditions. 
Prediction of a turbulent recirculating swirling flow using this code 
will be accomplished in Chapter VIII of this study. 


CHAPTER IV 


TURBULENCE MODELING OF TURBULENT FLOWS 
4.1 Introduction 

Most flows occuring in nature and in engineering applications are 
turbulent. Among those which an engineer commonly considers are: those 

which are unconfined by walls, such as the smoke plume from an indus- 
trial stack; flows near one wall, such as that in the vicinity of a 
turbine blade and an aircraft wings; and those which are confined duct 
flows. Those flows are characterized by hi§?>. levels of fluctuating 
vorticity. 

There are many books which approach turbulence in descriptive and 
experimental ways. Among the older texts is Townsend's book (53), which 
analyzes a wealth of experimental information with deep physical 
insight. Tenneks and Lumley (54) described a wide variety of fluid 
phenomena and the basic tools of turbulence theories, while Launder and 
Spalding (55) have provided in their book the main concepts of turbu- 
lence modeling and have shown what progress had been made and what 
problems remained. Finally, Wolfgang, and Bradshaw at al. (56, 57) 
provided two recent books which present a lecture series on prediction 
methods and turbulence modeling of turbulent flows. 

At the same time, many papers have given a considerable contribu- 
tion to turbulent flow understanding. Launder and Spalding (58) 
reviewed the problem of making numerical predictions of turbulent flow. 
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A turbulence model that permits predictions of both near wall and free 
shear flow was developed and evaluated by experimental data, while 
Saffman (59) eliminated some uncertainty regarding boundary conditions 
imposed at or near a solid boundary from turbulence model equations. 
Launder et al . (60) gave particular attention to the approximation of 
the pressure-strain correlations, and presented numerical solutions of 
the model equations for a selection of strained homogeneous flow and for 
two-dimensional inhomogeneous shear flows. Finally, a phenomenological 
theory of turbulence was given by Golowin (61). For further reading 
see, for example, References (62) through (70). 

4.2 Turbulent Flow Structure 

Turbulent flow of an incompressible fluid may be considered as a 
super position of eddies of different scales. The motion of large 
eddies which obtain energy from the mean flow is unstable and results in 
smaller eddies which in turn feed even smaller ones. In other words, 
the energy cascades from eddies of low frequency to those of high 
frequency. This process continues until the kinetic energy of turbu- 
lence dissipates (viscous dissipation). If there is no continuous 
external source of energy, the turbulence will decay. However, the 
interaction between the turbulent stresses and the velocity gradient in 
the mean flow prevents the turbulence from decaying by supplying the 
energy through the frequency spectrum. This process is termed 
'production 1 and the ratio of the production of energy to its viscous 
dissipation is a measure of the state of equilibrium of the flow. 

The eddy size at a given location in the flowfield is determined 
mainly by the flow domain and the flow local viscosity. For 
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sufficiently large Reynolds numbers, the small scale turbulence fluctua- 
tions are expected to be statistically isotropic (the statistical 
properties of the turbulence have no preference for direction) and to 
exhibit a universal scale similarity in their spectral behavior. How- 
ever, in the near wall region where the local Reynolds number decreases, 
the fine scale turbulence is no longer isotropic, regardless of how 
large the global Reynolds number becomes (71). The flow is not turbu- 
lent but viscous within the smallest eddies, and molecular effects are 
dominant (44). The gradient of the small scale eddies plays a dominant 
role in problems of heat and mass transfer by turbulence (72). 

4.3 Turbulence Modeling Approach 

From the discussion of the previous section, it can be concluded 
that the turbulent flow consists of large scale phenomena and small 
scale ones, and the characteristics of these scales are different. 
Therefore, it is very difficult to simulate them by the same approach. 
Large eddy simulation methods have been found useful in the study of 
large scale phenomence (73-76), and the small scale ones have been 
investigated by means of turbulent closure models (55, 77, 78). The 
former is a method midway between exact simulation and time-averaged 
modeled computation, which computes the large eddies explicitly and 
models only the small ones. Detailed descriptions of this method is 
beyond the scope of this study. The latter approach will be presented 
briefly in the next paragraphs. 

The most widely used turbulence closure models today are the eddy 
viscosity models and the Reynolds stress model. The first approach is 
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based on the eddy viscosity concept, which relates the Reynold stresses 
to the mean velocity gradient via: 


— i — r 3U. au . p 
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(4.1) 


where 6^ is the kronecker delta, k is the kinetic energy of tubulence 
(defined as 1 /2u ^ u ^ ) and y t is turbulent viscosity. 

The forms of this equation appropriate to an axisymetric flow are 
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From inspection of these equations, it is evident that for fully 
developed flow ( v and = 0 ) the normal stress components are equal, 
but in reality, the anisotropy is most pronounced. However, this defect 
is not too serious because the normal stress term is generally of minor 
importance in the Navier-Stokes equations for simple flows (e.g. thin 
shear layer flow). Also, the eddy viscosity concept breaks down in flow 
regions where the shear stresses and velocity gradient have opposite 
signs. Once again, this defect has little influence on the overall 
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predictions because regions with negative eddy viscosity (mathemati- 
cally) are in most cases small (71). Another observation is that zero 
velocity gradient provides zero shear stress, however, in flows where 
the asymmetry is pronounced, the positions of zero shear stresses and 
zero velocity gradient do not coincide (79). Generally, there is a 
limit to the applicability of models of this type, but the seriousness 
of their inherent defects will vary from one flow case to another. 

In the second approach (Reynolds stress model), individual turbu- 
lent stresses are represented by partial differential equations which 
may be derived from the Navier-Stokes equations. It should be noted 

that these equations contain correlations of one degree higher order 

i i _ i \ r 

correlations (e.g. u^u^ contains u.-UjU^ ) and other unknown functions of 
the velocity and pressure correlations as described in the next section. 
In order that these equations may be solved, all such correlations have 
to be modeled in terms of variables which are known or determinant, and 
the accuracy of the turbulence model is strongly related to this model- 
ing process. Thus a model of this type requires a large number of 
equations to be solved, which makes the model costly and difficult to 
use. 

Generally, the following steps are taken for the construction of a 
turbulence model (80): 

1. A set of properties that is sufficient to characterize the 
turbulent motion is chosen. 

2. A set of algebraic or differential equations which govern the 
distribution of these properties throughout the f'owfield is derived. 
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3. Constants or functions which are used to model some of the 
parameters are deduced from experimental data or from more fundamental 
theories. 

The univarsality of a turbulence model is controlled by the number 
of constants that need to be adjusted and by the number of differential 
equations that need to be solved. The greater the number of differen- 
tial equations and the lesser the number of constants that the model 
employs, the greater is the universality of the model. However, this 
increases the complexity of the computational problem and complicates 
the procedures of optimizing the set of constants. Complete universal- 
ity will never be achieved (80). 

4.4 A Brief Review of Turbulence Models 

A number of turbulence models have been proposed over the past few 
years. These models consists of either zero, one or multi-equations 
which when solved with the time-mean flow equations allow calculation of 
the turbulence properties and simulate the behavior of real fluids. 
Generally, most of the suggested models concentrate on the approximation 
of the Reynolds stress tensor which appears in the time-averaged Navier- 
Stokes equations. A brief description of selected zero, one, two and 
multi-equation turbulence models is presented in the following sub- 
sections. 

4.4.1 Prandtl Mixing Length Model 

Flows in a turbulent boundary layer are usually analyzed in terms 
of the Prandtl hypothesis, which relates the eddy viscosity directly to 
the local velocity gradient via: 
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The axisymetric form of this equation may be written as: 
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where 2 m is the mixing length which has to be prescribed empirically. 

The mixing length is usually expressed as a fraction of the width of the 
shear layer or mixing layer, and then generally assumed constant across 
the layer in regions remote from walls. The major shortcoming of this 
model is that it makes no allowance for diffusion and convective trans- 
port of turbulent viscosity. 


4.4.2 Energy Model [Prandtl 's 1945 Model] 


In order to eliminate the aforementioned shortcoming of the mixing 
length hypothesis, the eddy viscosity is expressed in terms of the 
length scale and the kinetic energy of turbulence k via: 

M t - C^pk l/2 Z (4.4) 

where C p is an empirical constant and i is a length scale of turbulence 
which is proportional to the size of the energy-containing eddies. 

Again, the length scale is a function of the flow geometry. The kinetic 
energy of turbulence is determined from the solution of a differential 
equation which expresses its transport. This model, with as a fixed 
empirical constant, is applicable only to flows or flow regions where 
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the local Reynolds numnber is sufficiently high (78). Hence, References 
(82-84) proposed low Reynolds number versions of the model by adding an 
exact viscous diffusion term to the k-equation and by replacing some of 
the model constants with functions of the turbulent Reynolds number. 
However, these new modifications have failed to predict strongly accele- 
rating boundary layers. Generally, this model works well only in rela- 
tive simple flows. 

4.4.3 Two-Equation Models of Turbulence 

These models employ the same formula for turbulent viscosity which 
was described in the previous model. However, in these models both the 
turbulence kinetic energy and the length scale are determined by solving 
their transport equations. Since the k field may be obtained from the 
solution of the k-equation, any combination of Z = k% n is possible to 
provide a mixing length field without solution of its transport equa- 
tion. Different combinations of m and n values give different models; 
examples of this: 

n = -2 and m = 1 gives k-u model 

n = -1 and m = 3/2 gives k-e model 

The latter gives the most widely adapted variable which is the dissipa- 
tion rate e {= k^/&). Both the dependence of the k-equation on e and 
the relative simplicity of deriving the e-equation have provided this 
parameter many advantages over the other combinations of n and m. 

Launder and Spalding (55) described the various two-equation models at 
length in terms of their derivations, applications and limitations. The 
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k-e model is the subject of the next chapter, however, discussion of the 
other two-equation turbulence models is beyond the scope of this study. 

4.4.4 Stress-Equation Models 


The models discussed so far lack the ability to account for the 
transport of the individual stresses. This is because they assume that 
the local state of turbulence can be characterized by one velocity scale 


u i 


and that the individual stresses can be related to this scale. 

In order to account properly for the transport of the individual 

stresses, models were developed which employ transport equations for 
— r — t 

u.u. . These models are often called second-order-closure schemes. 

* J 

The Reynolds stress transport equations for a fluid of uniform 
density p and with diminished buoyance effects may be expressed in this 
form: 
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or in words: 

Convection transport - 


■[diffusion transport.] 

• [stress production] + [pressure strain] 
■ viscous dissipation 


~ r Z 

Summation over the three equations for the normal stresses u^ provides 
the equation for the turbulent kinetic, energy k. The right hand side 
of Equation (4.5) contains a number of correlations of turbulent , 
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quantities. These correlations must be prescribed before the above 
equation is employed in a model of turbulence motion. Usually, the 
diffusion, pressure-strain and dissipation terms are approximated by 
models relating them to the Reynolds stresses, the mean velocity field 
and the length scale. Below are some approximation techniques which 
have been accepted by many workers for the turbulence quantities: 

1. For isotropic flow, the dissipation term can be approximated by 
2/3 e (Sjj where 5^- is as described in Equation (4.1) and e is the 
turbulence dissipation. 

2. The divergence of the fluctuating velocity transport equations 
derives the Poission's equation for pressure, which cart be used to 
eliminate the pressure-strain term (for example in Reference 60). 

3. References (60), (79) and (85) suggested different models for 
the diffusion term. 

Further details on the precise forms of the above turbulence 
quantities modeling is available in Reference (86). 

Taking into consideration the above approximations, a variety of 
models have been proposed (60, 87-89). A relatively simple model which 
was proposed by Launder and his co-workers (60, 89) is displayed here 
with the same previous assumptions. 
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where C$, and are constants that need to be prescribed. P^j is 

— r_r 

the stress production of component u.u. defined in Equation (4.5) and P 

' J 

is the stress production of the kinetic energy as defined in the 
transport equation of k. As mentioned in the previous section, these 
models consist of a large number of differential equations whose 
solution is costly. Hence the use of such models is limited. 

r 

4.4.5 Algebraic Stress Models 

An algebraic stress model is a simplified version of the Reynolds 
stress model. It has been derived by simplifying the differential 
transport equations such that they reduce to algebraic expressions which 
are easy to solve, but still retain most of their basic features. The 
differential equations are converted into algebraic expressions by elim- 
inating the gradients of the dependent variables which appear in the 
transport equations. The resulting general equation which corresponds 
to Equation (4.5) is: 

u i u j = P+e (C^l) % 6 ij {C 2 P + e(C l' i)} + P ij ( 1_C 2 ^ ( 4 -') 

where P and P^- are as defined in Equation (4.6). The background to two 
equations is described at length elsewhere (86, 90). Application of a 
variety of these models has been reported in References (78), (91), and 

(92). 

4.5 Summary 

The present chapter provides a brief description of the structure 
of turbulent flows and the mathematical models which are expected to 
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represent the flow behavior. The first three sections consist of a 
survey of the contributions of many workers to the understanding of the 
characteristics of the turbulent flows, and with the turbulent flow 
structure and its modeling approaches. Various turbulence models are 
discussed in the fourth section. The Prandtl mixing length model was 
presented to demonstrate the non-equation models, and the energy model 
represents the one-equation models. The k-e turbulence model is among 
the two-equation models which are recommended for many flows. As the 
flow becomes more complex, a higher order model is required, thus the 
Reynolds stresses model and algebraic stress model were presented. 
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CHAPTER V 


THE k-e TURBULENCE MODEL AND RECENT EXTENSIONS 


5.1 The k-e Turbulence Model 


The k-e model has been used by many workers (2, 43, 55, 93-95) to 
provide realistic predictions of the velocity fields and the turbulence 
parameters. This model involves two transport equations for turbulence 
characteristics. One of these equations governs the distribution 
throughout the field of the turbulence energy k, which measures the 
local kinetic energy of the fluctuating motion; the other governs the 
turbulence energy dissipation rate e. The k equation can be obtained 


from the exact, transport equation of the Reynolds stresses u. u. by 

* 0 

summation over the three equations for the normal stresses. The result- 
ing equations for the axisymmetric cylindrical flows are as shown below: 

u. 


9 k 9 k _ 1 9 r ^ e 9 k -| 9 | - M e 9k \ ~ 

p u 87 + pv a? - 7 3? [r ^ 37 ] + aV * G ‘ pe 


9e , 9 £ _ 1 9 r 1 e 9s i ,9 r^e 9e-i . r r n, 

p u a* * pv 3r " 7 57 [r T s? J + li S ] + C 1 G e/k 

£ £ 


(5.1) 


- C2 p e/k 


The values usually adapted for the empirical constants C^, C 2 , cr^ and a E 
are those quoted in Reference 55; they are listed in Table IV. The 
generation term is: 
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where w e is the sum of turbulent and laminar viscosities. Knowledge of 
k and e allows the length scale z and turbulent viscosity uj. to be 
determined from which the turbulent shear stresses can be evaluated via: 


w t * c u p k2 / £ (5.2) 

and Equation 4.2. 

To account for the wall effect in the near-wall regions, equations 
are introduced to link velocities, k and e on the wall to those in the 
near-wall region. These equations, called "wall functions", are intro- 
duced and used in finite difference calculations at near-wall points. 
They occur in the momentum equations and k-generation terms, and their 
implementation is discussed elsewhere (55) together with appropriate 
near-wall e specification. The main assumptions considered in these 
formulations are: 

1. Convection and diffusion of turbulence kinetic energy are 
nearly negligible. 

2. Eddy viscosity is isotropic. 

3. The sublayer extends up to a value of wall function parameter 
y + equal to 11.63. 

4. The length scale is proportional to the normal distance from 
the wall . 

5.2 The k-e Empirical Constants 

References (79) and (96) discussed in some detail the basis for 
choosing the k-e constants. Therefore, it may suffice to mention here 
that Cy is fixed by the requirement that in a constant stress layer 
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production and dissipation of e are balanced [x w /(pk) = Cg is 

determined by reference to the decay of grid turbulence, and Cj is 
chosen so that the von Karman constant equals 0.42, The 'Schmidt' 
numbers and a E are fixed by computer optimization in Reference 
(97). The main assumptions considered in determining the k-e constants 
are: the flow is not far from equilibrium so that the convection and 

diffusion are small compared to production and dissipation terms; and 
the turbulent Reynolds number is high. These assumptions have led to 
the nonuniversality of the model constants; even in certain fairly 
simple flows, some constants require different values. Usually, the 
extension of the model applicability is achieved by algebraic 
formulation of Cg and C p to account for a specified flow case. These 
formulations will be discussed later. 

5.3 Shortcomings of the k-e Model 

The width of applicability of this model was tested earlier by 
Reference (59), who considered nine substantially different kinds of 
turbulent flow. The model demonstrated the capability of prediction of 
both near-wall and free shear flow phenomena without adjustments to 
constants or functions. However, in some flow cases this model has some 
shortcomings, a brief summary of which is provided below. 

Habib and Whitelaw (35) found that the velocity minimum and maximum 
on the centerline for confined swirling flows with coaxial jets are not 
correctly predicted. However, the measured values of the two components 
of the rms value of the velocity fluctuation and the kinetic energy are 
in excellent agreement with the predicted values. This led to their 
support of Gibson (91) that the assumption of a constant value of C u is 
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inappropriate in curved flows and they remark, that empirical adjust- 
ments to constants may produce better results. 

Syed and Sturgess (98) predicted the rms value of axial velocity 
fluctuation, turbulence parameters, velocity profiles and turbulence 
shear stresses for coaxial jets in confined sudden expansion flows. A 
maximum error of about 27 percent was obtained in the rms values. The 
maximum discrepancy in the predicted velocity profiles was found at the 
centerline. The maximum level of turbulent shear stress was predicted 
as occurring at a lesser radius than was measured. This led to the 
conclusion that the dissipation rate is not represented correctly and/or 
the convection of shear stress is not negligible and its value is not 
based mainly on the local properties. 

Measured and predicted streamline contours, velocity and turbulence 
intensity profiles for isothermal flow in axisymmetric models of 
combustor geometries were compared by Whitelaw and Green (99). It was 
seen that the turbulence modeling errors are prevalent in regions of 
recirculation, flow development, and in the near-wall region. 

The finding of Sturgess and Syed (45) for isothermal flow in widely 
spaced coaxial jet, diffusion-flame combustors was that the transport of 
turbulent stress is involved in producing the local centerline axial 
velocity peak. This peak cannot be predicted by the eddy-viscosity 
model because the turbulence characteristics are derived from the 
velocity profiles and eddy viscosity. Such models are incapable of 
calculating the correct flow characteristics in regions of negligible 
velocity gradient, particularly if transport of turbulent shear stress 
of a sign opposite to that appropriate for the local velocity gradient 
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is significant. Therefore, Reference (45) recommended to take into 
account the transport of turbulent shear stress. 

5.4 Previous Proposed Extensions 


The above mentioned problems have commanded considerable attention 
from many investigators. Most of the investigators have agreed that 
some of the model shortcomings can be eliminated by algebraic formula- 
tion of C 2 and C y . Usually, these corrections are expected either to 
account for low Reynolds number or streamline curvature effects, or 
both. 


Early, Jones and Launder (97) added a new term. 


7 ,3k 1/Z i2 

- 2 “ ( iy ) 


to the k-equation to represent the energy dissipation rate at the wall. 
Also, another term, 

2 v ( — o) 

ay 

was added to the ^-equation in order for the peak level of the turbulent 
kinetic energy to agree with experimental data. Both C 2 and were 
modeled to account for low Reynolds number effect via: 


C y = 0.09 exp [-2. 5/(1 + 1^/50)] 
C 2 = 2.0 [1.0 - 0.3 exp (-R^)J 


(5.3) 


where R^ = p k 2 /(eu ). 

This model was tested by Jones and Launder (100) as applied to 
isothermal low Reynolds number pipe flow, and wall boundary layer with 
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streamwise pressure gradient and wall injection. The model predictions 
were not completely in accord with the experimental data. 

Launder et al. (101) incorporated the streamline curvature effects 
on the turbulence structure in the k-e model. The direct effect of 
curvature in the model is limited to a single empirical coefficient 
whose magnitude is directly proportional to a Richardson number based on 
a time scale of the energy-containing eddies. For boundary layer flow 
types, the proposed C 2 formulation is: 

C 2 = C 2 (1 - 0.2 R.) (5.4) 

I 

where C 2 is the value before the correction, and the Richardson 
number is defined as: 

R-j ~ ^2 ^2 yp (fw) (5.5) 

e r 

where w is the swirl velocity. The model was applied to a variety of 
boundary layer flows developing over curved and spinning surfaces; 
satisfactory agreement with experiments was obtained. On the other 
hand. Reference (102) solved Reynolds equations using the k-e model in 
the separated flow region, and both the k-e model and the k-e-u ' v T 
models in the region downstream of the separation. He showed that in 
the redevelopment region the three-equation turbulence model does not 
offer worthwhile advantages over the two-equation k-e model. 

Ljuboza and Rodi (103) replaced the empirical constant C w by a 
function which was derived by reducing a model from the Reynolds stress 
transport equations to algebraic expression by introducing simplifying 
assumptions about the convection and diffusion terms. Good results were 
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obtained when this model was applied to a wall jet in stagnant surround- 
ings, or in a moving stream. 

Hanjalic and Launder (104) introduced into the e-equation a term 
proportional to 

8u.j 3Uj 

k 3xJ 3x m e i jk e lmk 

where e-jj^ stands for the third-order alternating tensor. For thin 

9 u 2 

shear flows this term reduced to k (-r-) . This correction is expected 

°y 

to promote a higher rate of dissipation for irrotational strains as 
opposed to for rotational strains. Indeed, they found that this minor 
modification brought substantial improvement to the prediction of 
boundary layers in an adverse pressure gradient. 

The streamline curvature correction proposed by Reference (77) was 
derived from the Reynol ds-stress equations, as described by Reference 
(91). Also, they adapted a modification which accounts for the prefer- 
ential influence of normal stresses and the dissipation of turbulent 
energy; this modification was originally proposed by Reference (105) for 
strong accelerated flows. This correction has led to greatly improved 
agreement with experimental data when applied to the recirculating flow 
in an annulus, and a plane twin-parallel jet in still air. 

Karasu (5) used the turbulent viscosity u e , which is calculated 
from the k-e model, in the solution of the axial and radial moments of 
swirling pipe flow. The tangential eddy viscosity y r x was employed in 
the solution of the tangential momentum; the expression used is given by 
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(5.6) 
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where k is a constant considered to be 0,042 and \i is the laminar 
viscosity. His prediction for pipe swirling flow that exhibits combined 
vortex were generally comparable with the experimental data. However, 
the prediction of the axial velocity at the initial location of pipe 
flow that exhibits so'iM body rotation was poorly compared with the 
measurements. 

Finally, good agreement with experimental data was obtained (106) 
when a modified k~e model was applied to fully developed pipe flow. The 
modified version can be used without introducing the empirical wall 
function formulas. 


5.5 Closure 

A brief description of shortcomings and recommended corrections to 
the standard two-equation k-e turbulence model suggested by previous 
investigators are presented. It has been found that the model predic- 
tive capability can be increased by a simple empirical modification of 
the model governing equations and this modification works only well for 
the flow under consideration. Thus the universality of the model has 
never been obtained. 
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CHAPTER VI 


TURBULENCE PARAMETERS DEDUCTION FROM SWIRLING 
CONFINED FLOW MEASUREMENTS 

6.1 Introduction 

Specification of the turbulent stresses is required in order to 
solve the Reynolds equations for time-mean velocities. This is usually 
obtained via the turbulent viscosity y t . Chapter V exhibits some simple 
models for the specification of u t , especially in boundary layer flows, 
however, these models are not useful for recirculating flows. 

Currently, two-equation models are popular in these circumstances. The 
present work emphasizes the k-e version of the two-equation models in 
the context of swirling recirculating flows. As mentioned earlier, this 
model solves differential equations for the turbulence energy, k, and 
the turbulence dissipation, e. From these the turbulent viscosity can 
be specified. 

From the above introduction, it is foreseen that the precision of 
the predictions strongly depends on the accuracy of the turbulent 
viscosity values. The evaluation of a turbulence model and hence the 
values can be performed through the use of analytical inverse methods 
which are described in the next section. 
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6.2 Analytic Inverse Methods 


Evaluation of the accuracy of a turbulence model may be obtained by 
comparing the calculated time-mean values and turbulence parameters with 
experimental data. Alternatively, Equation (4.2) may be used to make 
deductions directly from experimental data. Velocity gradients are cal- 
culated from experimental time-mean velocities using a suitable finite 
difference technique. Once the velocity gradients are obtained the 
turbulent viscosity components can be calculated via Equation (4.2). 
Implementing the Prandtl mixing length model, the mixing length can 
be deduced from Equation (4.3) where Equation (4.2) provides the value 
of u rx . The other length scale calculation follows from replacing by 
u rx in Equation (4.4) and setting C u to its conventional value of 0.09. 
Resultim length scale values are presented in Section 6.3. 

In using Equations (4.2) and (4.3), lack of detailed axial time- 
mean velocity gradients has led to these terms being dropped. In non- 
swirling and moderately swirling flows, they are usually small, except 
near recirculation zones. At high swirl, the central recirculation zone 
develops and this simplification cannot faithfully be invoked (107). 
However, it will be only in quite small regions of the flow that errors 
of any significant magnitude might occur. 

6.3 Findings and Discussion 


The present chapter was initiated to analyze the fluctuating and 


time-mean turbulence data presented in Reference (10). The data are for 
nonswirling and swirling flows in an axisymmetric test section with 
expansion ratio D/d = 2, which may be equipped with a contraction nozzle 
of area ratio 4 located at L/D = 2. The swirl is imparted to the flow 
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by means of an adjustable-angle vane swirl er (3). The Boussinesq's eddy 
viscosity was used, as explained in the previous section, to calculate 
turbulent viscosity component from each measured Reynolds stress and its 
corresponding time-mean velocity gradient term. Also, from this mean 
component of turbulent viscosity the length scale of energy containing 
eddies was obtained using the measured kinetic energy of turbulence with 
the parameter given its standard value, 

6,3,1 Assumptions and Parameters Normalizations 

Based on the available data, curves are usually fitted to the 
spatial distribution of measured time-mean values before further analyt- 
ical manipulation. However, in case of complex flows, it is very 
difficult, or it may be impossible, to find a general correlation of the 
data along the flowfield. Therefore, the data were hand smoothed and 
then the values are read in from the corrected curves. The data consid- 
ered were for swirl vane angles of 0, 45 and 70 degrees. The radial 
velocity gradients were determined analytically from the corrected data. 
The axial velocity gradients were ignored because their values were 
unavailable and are expected to be considerably less than the radial 
gradient values in most flow situations studied here. 

The turbulent viscosity was normalized with respect to the swirler 
inlet uniform axial velocity u 0 (deduced independently from a measure- 
ment upstream of the swirler), fluid density and combustor diameter. 

The turbulence kinetic energy was normalized with respect to u 0 2 . The 
length scale and mixing length were nondimensionalized with respect to 
the combustor diameter. The normalized values of the viscosity in the 
r-x plane u rx , the length scale and the kinetic energy are displayed in 
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Tables V through X. Also, the other viscosity components are tabulated 
as a fraction of u rx values (a^ = y rx /u-jj) and displayed among these 
tables. The tabluated y rx and length scale were computer plotted and 
are shown in Figures 5 through 10, 

6.3.2 Effects of Swirl on Turbulence Parameters 

The swirl vane angles 4> considered in this study are 0 (swirl er 
removed), 45, and 70 degrees. Part (a) of Figures 5 through 7 displays 
the nondimensionalized turbulent viscosity u rx and part (b) shows the 
normalized length scale z as defined in Equation (4.4) for <t> = 0, 45, 
and 70 degrees, respectively. Also, the corresponding values of these 
parameters are displayed in parts (a) and (b) of Tables V through X, 
respectively. It should be noted that the allowed maximum range of the 
values plotted in Figures 5 through 7 was 2, and the data beyond this 
range was assigned the range maximum values. The actual values of the 
affected points can be found in Tables V through VII. Parts (c) through 
(f) of these tables exhibit the other turbulence parameters such as 
turbulence energy k and other turbulent viscosity components. 

The y rx profiles shown in part (a) of Figures 5 through 7 reveal 
that the y rx values marginally increase as the swirl strength increases 
from zero to a moderate value and sharply increase beyond this range. 
Comparison of these figures with the mean velocity plots presented in 
Reference (10) indicates that the peak of u rx may or may not coincide 
with the peak of the velocity gradients. Inspection of part (b) of 
Figures 5 through 7 shows that the normalized length scale profiles are 
almost similar to the corresponding u rx profiles for each swirl case and 
therefore, the swirl dependency of length scale is revealed. Parts (a) 
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and (b) of Tables V through VII tabulated the plotted values of the 
normalized turbulence viscosity and the length scale, respectively. 

Part (c) of Tables V through VII displays the mixing length values 
obtained from Equation (4.3), with the axial velocity gradients ignored, 
using the curve fitted data of Reference (10). Comparison between these 
tables reveals that the mixing length values are slightly higher for 
swirl vane angle <f> for 45 degrees than for $ = 0 or 70 degrees. How- 
ever, the presented values for <p = 0 degrees are smaller than the corre- 
sponding ones for 70 degrees. It is expected that these findings are 
due to the effect of swirling motion of the flow on the eddy size 
distribution. 

Part (d) of Tables V through VII contrasts the kinetic energy of 
turbulence distributions as the swirl strength increases. The tables 
show that results for moderate swirl are slightly higher than the corre- 
sponding ones for nonswirling flow. The most dramatic effect of swirl 
is seen at high swirl values. The turbulence energy maximum values 
occur in regions of recirculation and regions of high stress. 

The degree of nonisotropy has increased in magnitude, consistent 
with the increase in swirl strength throughout the entire flowfield, as 
can be seen from part (e) of Tables V through VII. Also, it can be 
noticed that the turbulent viscosity in the rx-plane may be higher or 
lower than its corresponding value in the other planes depending on the 
location of the point in the flowfield. 

Generally, the effect of increase of swirl strength is the consid- 
erable increase in all the turbulence parameters considered--that is, 
increase in rx-viscosity, kinetic energy of turbulence, length scales, 
and degrees of nonisotropy. The peak location of these parameters 
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depends on the parameter under consideration. Overall, the swirl 
produces a large-scale effect on the turbulence parameters and, there- 
fore, its effects should be considered in the turbulence modeling. 

6.3.3 Effects of Strong Contraction Nozzle 
on Tubulence Parameters 

Comparison of results of the previous sub-section with the corre- 
sponding ones of Abujelala et al. (108) shows that the general results 
which were obtained from the corrected curves are similar to those 
obtained directly from the scattered data. Therefore, a decision was 
made to employ the measured values in the calculations involved in this 
sub-section. The normalized turbulent viscosity in the rx-plane and the 
length scale for swirl vane angles of 0, 45 and 70 degrees, and with a 
strong contraction nozzle of area ratio 4 located at L/D = 2 are 
presented in parts (a) and (b) of Figures 8 through 10, respecti vely . 

The nozzle diameter is half the combustor diameter. Some of the plotted 
points were assigned the maximum values allowed if the normalized values 
were larger than 2. Actual values are included in the corresponding 
tables. The effect of the contraction nozzle on turbulence parameters 
was tabulated and is presented in Tables VIII through X. The construc- 
tion of these tables is similar to Tables V through VII, 

Comparison between Figures 8 through 10 and Figures 5 through 7 
reveals that the nozzle tends to increase the u rx and z values in the 
shear layers and to decrease them elsewhere. It is suggested that this 
is because of the high velocity gradients which are generated by accel- 
erating the flow near the centerline. The similarity between the pro- 
files shown in parts (a) and (b) of Figures 8 through 10 was unaffected 
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by the presence of the nozzle. The data of these figures were obtained 
from parts (a) and (b) of Tables VIII through X. 

Parts (c) through (g) of Tables VIII through X show that the con- 
traction nozzle generally changes the turbulence characteristics. This 
effect increases as the swirl strength increases. Parts (c) and (d) of 
these tables reveal that the nozzle tends to increase the mixing length 
and the kinetic energy in the region of recirculation, and to reduce 
their values near the wall. This is because the presence of a contrac- 
tion nozzle accelerates the flow near the centerline and decelerates it 
close to the wall. Parts (e), (f) and (g) of Tables VIII through X show 
a high degree of nonisotropy due to the presence of the nozzle. This is 
expected because of the streamline curvature effects which increase as 
the contraction increases. Summarizing, the presence of a strong con- 
traction nozzle change the turbulence characteristics dramatically, and 
its presence should be taken into consideration in turbulence modeling. 

6.4 Closure 


The data of confined swirling flow obtained from a hot-wire using 
the six-orientation technique were numerically analyzed. The effects of 
swirl strength and a strong contraction nozzle on the turbulence param- 
eters were evaluated. It was found that with the swirl strength and the 
contraction have strong effects on the turbulence parameters. 

Generally, the most dramatic effect of tlie increase of swirl is the 
considerable increase in all the parameters considered. The presence of 
a strong contraction nozzle tends to increase the parameter values in 
regions of acceleration where large radial velocity gradients occur and 
to reduce them in the deceleration region near the outer boundary. 
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CHAPTER VII 


EXTENSION OF THE k-e TO CONFINED SWIRLING 
RECIRCULATING FLOWS 


7.1 Introduction 


Turbulent swirling recirculating flows are used in many practical 
situations. They provide the fundamentals of physical processes occurr- 
ing in aircraft combustors and in 'industrial furnaces. These types of 
flow are characterized by strong streamline curvature, complex eddy 
structures and high turbulence intensities. The practical need for 
computation of this category of flows has encouraged many workers to 
attempt to develop turbulence models that r;«n provide accurate informa- 
tion about such flows. However, most of these attempts have provided 
models that are only applicable for a certain flow, and these models 
must be modified to accommodate flows of greater complexity. 

The different models may conventionally be distinguished by the 
number of turbulence properties which appear as dependent variables of 
the differential equations, as described in Chapter IV. The most widely 
used model is the k-e turbulence model which was briefly described in 
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recirculating confined swirling flow (2): it provides inaccurate predic- 
tions for the size and strength of the recirculation zones. It is 
expected that the turbulence model over-predicts the dissipation of the 
kinetic energy of the flow, therefore, the centrifugal forces produce a 
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shorter central recirculation core than the measurements show. Hence, 
in order to extend the k-e model to confined swirling flows, some 
corrections concerning the e-equation have to be made. The general 
approach, as shown in Chapter V, is to consider the streamline curvature 
effects through the modeling of the constant C 2 . This follows in the 
next section. 


7.2 Previously Proposed Extension and 
Its Invalidity 


Reynolds (109) defined the Richardson number for swirling flows as: 


R. = 


2 (r "> 

r 

(M) 2 + [ r — (-)3 2 

r' L 3r v r' J 


(7.1) 


An extensive preliminary study was conducted to evaluate this equation 
and Equation (5.5) as applied to turbulent swirling recirculating flows. 
The distribution of predicted C 2 and mean velocity values using the two 
expressions for the Richardson number along with Equation (5.4) were 
compared and the results show that Equation (7.1) provides better 
predictions than Equation (5.5). Therefore Equation (7.1) along with 
Equation (5.4) was selected for the evaluation of the inclusion of the 
streamline curvature effects through modeling C 2 . 

Equation (5.4) with the Richardson number as defined by Equation 
(7.1) was added to the e-equatlon of the k-e, and the flowfield for 
swirl vane angles of 45 and 70 degrees was predicted using STARPIC 
computer code which is described in Chapter IV. The value of C 2 was not 
allowed to be less than 0.1 and not greater than 2.4, with C y and the 
other model constants assigned the conventional values. The distribu- 
tion of C 2 values over the flowfield is displayed in Tables XI and 
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XII. Table XI displays a moderate swirl strength flowfield ($ - 45 
degrees), and Table XII exhibits the flowfield for a stronger swirl 
strength (4 = 70 degrees); the values not included in the above 
mentioned range have been replaced by a {-) sign. The corresponding 
velocity profiles for $ = 45 and 70 degrees are displayed in Figures 11 
and 1 2, respectively. Inspection of the predicted data associated with 
these figures reveals that the corrected e-equation under-predicted the 
dissipation of the flow energy. The dominant centrifugal forces drive 
the flow off the centerline. Similar results were obtained for the <j> = 
70 degrees. 

7.3 Parameter Optimization From Time-Mean Data 

It is clear from the previous section that the variation of the Cg 
parameter given in Equation (5.4) does not enhance the predictions of 
turbulent swirling recirculating flows, and it was decided to optimize 
the k-e model constants for the flow under consideration. It was 
mentioned in Section 5.2 that only and a e were fixed by computer 
optimization. However-* because of lack of the necessary experimental 
data, C p arid C 2 were also optimized in this study. The constant was 
dropped from the optimization because of the feeling that the k-equation 
produces sufficient accuracy in determining the k-values. An effective 
optimization code is the Hill Algorithm. Kuester (110) described this 
method at length, in terms of its mathematical formulation and practical 
applications. It need be described only briefly here. Therefore, it is 
described only briefly in this section. 

The optimization procedure is based on the direct search method 
proposed by References (111) and (112). This method has proven 
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effective in solving some problems where the variables are constrained. 
The procedure requi res a starting point that satisfies the constraints 
and does not lie in the boundary zones. The effect of a small stepping 
in the first variable on the objective function is evaluated. If an 
improvement in the objective function were obtained without violating 
the boundary zones or constraints, the step size is reduced, and the 
direction of movement is reversed. The next variable is in turn stepped 
by a small distance parallel to the axis. The same acceleration or 
deceleration and reversal procedure is followed for all variables in 
consecutive repetitive sequences until the convergence criteria is 
satisfied. 

The parameters C y , Cg, and a £ were optimized by minimizing the 
maximum absolute discrepancy between predicted axial and swirl velocity 
and the corresponding experimental values. These experimental values 
were obtained using a five-hole pitot probe (9). The optimum constants 
were obtained by an iterative computer procedure; the velocity fields, 
after the iteration limit was exceeded, were compared with the experi- 
ments. Because of the rather large computing time that would have been 
required if the iteration was continued until the convergence criterion 
was fully achieved, only 80 iterations were performed during each opti- 
mization cycle, with current flowfield values taken as initial values 
for the next optimization cycle in each case. The k-e constants optimum 
values obtained after 20 optimization cycles were found to be only 
marginally dep ndent on swirl strength. The values obtained for this 
wide range of swirl strengths (<f> = 45 to 70 degrees) are displayed in 
Table XIII. It should be noted that these values are not as good as the 
Table IV values in the prediction of nonswirling flows. 
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7,4 C y Formulation From Detailed Turbulence Data 

The analysis of the detailed turbulence measurement of (10) pre- 
sented in Chapter VI are being employed here to provide the background 
for producing an empirical formulation for the k-e C p parameter. 
Insepctlon of the normalized rx-viscosity „ rx and length scale i, shown 
in Figures 5 through 10 of that chapter reveals a strong similarity 
between the radial profiles. That is 

* l t / ( pu o D) 

is almost proportional to z/D. Since 

u t = c „ pkl/2 * 

or equivalently 

u t /(pu 0 D) = (C M k 1/2 /u 0 )(t/0) 

it follows that 

Ck^/u 

is approximately constant throughout the flowfield. That is: 

c v = Au o /kI/2 (7.2) 

Where A is a constant Independent of spatial position in the flow. 

Using an optimization procedure similar to that described in Section 

7.3, optimal values of A, Cg and o e have been deduced for the present 
situation: It was found that 
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A - 0.0083 


C 2 = 1.804 (7.3) 

o e = 1.455 

The flowfieids for inlet swirl vane angles of $ * 45 and 70 degrees 
have been predicted using the C p variation as given by Equation (7.2), 
with the optimized values of Equation (7.3) used in this and other 
equations occurring in the simulation. Parts (a) and (b) of Table XIV 
exhibit the values of predicted k and C p obtained in this manner for 
corresponding inlet swirl are angles of <J> = 45 and 70 degrees, respec- 
tively. Values vary throughout the flow domain. The predictions using 
the aforementioned C g variation are validated in the following sections. 

7.5 Predictive Capability of C u Formulation vs. 

Optimized Parameters 

The time-mean flowfields corresponding to the swirl vane angle 
<j> - 45 and 70 degrees, which were obtained by two different approaches, 
are compared in this section. Parts (a) of Tables XV and XVI display 
the predicted flowfield using the C y formulation for ij> = 45 and 70 
degrees, while parts (b) of these tables exhibit the predicted flowfield 
using the optimum k-E model constants for the same above swirl 
strengths. 

Generally, the comparison between Tables XV and XVI reveals that 
both methods provide very close predictions. However, the predicted 
centerline velocity values of parts (a) of these tables reduced 
considerably as compared with the parts (b) values, .and became close to 
the experimental values as will be seen in the next section. Thus the 
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C p formulation predictions are superior over the ones of the optimum 
values in the region of recirculation. 

7.6 Validity of Standard k-e Model and Formulation 
7.6.1 Time-Mean Axial and Swirl Velocities 

The predicted velocity profiles for swirl vane angles of 45 and 70 
degrees using the C M formulation of Equation (7.2) are compared with the 
five-hole probe (9) and hot-wire (10) measured data, as can be seen in 
Figures 13 and 14. respectively. Inspection of these figures reveals 
that very good agreement between the predicted results and the experi- 
mental data has been achieved. Note that the measured inlet profiles 
are plotted at the x/D = 0 location; in fact these are actually the 
values measured immediately downstream of the swirler, at the location 
x/D - -0.11. It is convenient to retain these values on the profile 
plots in the prediction study, although clearly results at the x/D = 0 
location are then not only directly comparable with the inlet station 
data of References (9 and 10), which are taken precisely at x/D = 0. 
These comments apply also to other plots given in the thesis. The 
discrepancies in the centerline velocity values at the peak of the 
central recirculation core are expected due to the assumption of 
a xi symmetry at that location. 

Figures 15 and 16 exhibit predictions similar to Figures 13 and 14, 
except the standard k-e model was used in these predictions. Inspection 
cf these figures reveals that inaccurate predictions for the size and 
strength of the recirculation zone were obtained by using the standard 
k-e model. The inaccuracy of the standard model can be easily seen by 
comparing the measurements of Yoon and Lilley (9), and Jackson and 
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Li 1 1 ey (19) for swirl vane angles of 45 and 70 degrees as shown in 
Figures 13 and 14 against the corresponding predicted velocity profiles 
displayed in Figures 15 and 16, respectively. However, as mentioned in 
the previous paragraph, the predicted results using the C p formulation 
has shown very good agreement with experimental data. Thus it is clear 
that the C M formulation has provided superior results over the standard 
k-e model . 

7,6,?. Zero Velocity Envelope and Time-Mean 
Centerline Velocity 

Figures 17 and 18 show a comparison of the predicted central recir- 
culation zone boundary using the standard k-e model and its modified 
version (variable C p ) with corresponding data (9, 10) for swirl vane 
angles of 45 and 70 degrees. The predicted central recirculation zone 
envelope using C y variation is in very good agreement with the data, 
with a small discrepancy near the peak of the recirculation zone core. 
However, using the standard model produces a less realistic shorter 
central recirculation zone. 

Predicted centerline velocity longitudinal variation using the 
aforementioned methods as compared with the measurements of References 
(9) and (10) is shown in Figures 19 and 20 for the same previous swirl 
strengths. Notice that the negative values are slightly over-predicted 
in absolute magnitude--a trend which was also found by others (35, 98) 
as reported in Chapter V. 

7.6.3 Turbulence Parameters 

The predicted kinetic energy, k and its dissipation rate, e using 
the standard and the modified version of the k-e model are compared with 
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the corresponding hot-wire data for swirl vane angles 45 and 70 degrees 
as shown in Figures 21 and 22, respectively. Both versions of the k-e 
model have provided qualitatively good predictions. The modified model 
shows closer predicted values of k and e to the experimental data 
(especially within the recirculation zone and far downstream of the 
inlet domain, respectively) than the standard model. 

7,6,4 Summary 

In summary, the C y variation has provided better predictions of the 
time-mean and the turbulence parameters than the standard model. Also, 
the size of the recirculation zone, especially for swirling flow and the 
centerline velocity are quantitively very well predicted by using the 
proposed formulation as described in Section 7.4 of this chapter. 

7.7 Validation of the Proposed C„ Formulation 
for Nonswirling and Pipe Flow Applications 

It was mentioned earlier. Section 7.4, that the C y formula was 
obtained via a simple optimization procedure. The optimum values were 
obtained by minimizing the difference between the predicted and measured 
axial and swirl velocity values for the flowfield of swirl vane angle = 
45 degrees. The validity of the proposed formulation as applied to a 
higher swirl strength flow field (70 degrees) has been shown in the 
previous section. Two more test cases are being Conducted in this 
section to show the generality of this model. A nonswirling flow 
through a circular duct with expansion ratio d/D = 2 and a swirling 
flowfield with d/D = 1 (pipe flow) are considered for this test. 

Figures 23 through 25 give the comparison between the predicted 
time-mean velocity fields and the corresponding measurements. Predicted 
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axial velocity profiles for nonswirling flow are compared with the 
experimental data of Yoon and Lilley (9), as seen in Figure 23, Figures 
24 and 25 exhibit the comparison between the predicted axial and swirl 
velocity profiles and the measurements of Jackson and Lilley (10) for 
swirl vane angles 45 and 70 degree fields, respectively. An examination 
of these figures reveals that the predictions are in very good accord 
with experimental data, particularly when the experimental scatter is 
borne in mind. This finding along with the previous section conclusion 
lead to the generality of the proposed C M formula for applications to a 
circular duct flow with and without abrupt expansion, 

7.8 Closure 

An evaluation of the standard k-e model and its previously- 
available extensions as applied to turbulent reci rculating confined 
swirling flow has been presented. Neither the standard model nor its 
extensions show effective capability for predicting turbulent recircu- 
lating swirling flows. Recent experimental data on swirling confined 
flows, obtained with a five-hole pitot probe and a six-orientation hot- 
wire probe, are used to obtain optimum values of the turbulence param- 
eters , Cg and and a special C y formulation for swirling flows', 
respecti vely. General predictions of moderately and strongly swirling 
flows with these techniques are more accurate than the predictions with 
the standard or previous simple extensions of the k-e turbulence 
model. The predicted centerline velocity using the C p formulation has 
been found superior to the correspondi ng velocity given by the optimum 
values. The application of the formula to nonswirling flow and pipe 
flow has shown the generality of this model. 
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CHAPTER VIII 


GENERAL PREDICTIONS FOR CONFINED TURBULENT 
SWIRLING FLOW 

8.1 Scope and Method of Approach 

This chapter first discusses predictions which deal with the flow 
through the inlet swirler with vane angles set to <j> = 45 and 70 
degrees. Various types of inlet profiles assumptions are considered and 
the similarity and differences in the ensuing flowfield predictions are 
noted. Predictions are then exhibited for a range of swirl strengths 
correspond to i> ~ 0, 45 and 70 degrees using measured inlet axial, 
radial and swirl velocity profiles in each case. Downstream nozzle 
contraction ratio and inlet domain expansion ratio effects are inves- 
tigated. The ensuing flowfields are characterized via velocity profiles 
and streamline patterns, and illustrate the large scale effects of inlet 
swirl, expansion ratio and outlet nozzle on flowfields. 

The prediction procedure starts with the solution of partial 
dif ferential equations of conservation of mass, momentum in x, r, and 0 
directions, turbulent kinetic energy and its dissipation rate, which 
govern two-dimensional axisymmetric steady flow via the STARPIC computer 
program, the solution techniques which have been implemented in this 
computer program are described at length in Chapter IV. All the results 
to be given in Sections 2 through 6 of the present chapter are obtained 
via the formulation model as described in the previous chapter, see 
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Section 7.4. Similar predictions using the standard k-e model are pre- 
sented for flowfields with D/d = 2, 1.5 and 1. The former was obtained 
earlier (2} and shows the same general trends for the effects of swirl 
and downstream contraction nozzle on the flowfield, it is not shown 
here. The predictions of the other two flowfields are discussed in 
Section 8.5. This chapter emphasizes the more advantageous 
formulation model, which gives a more accurate simulation of the 
flowfield. 

The code is operated in the manner described in Chapter IV, with 

the desired parameters being computed to convergence with a 23 x 21 

nonuniform grid covering a flow domain of length 4 times the combustor 
chamber diameter for sequential swirl vane angles of 0, 45, and 70 
degrees. 

The open ended flowfield with expansion ratios of 1.5 and 1 are 
covered by 23 x 15 and 23 x 9 nonuniform grid systems, respectively. It 
should be noted that the refinement of the grid was necessary only for 
highly and moderately expanding flow. The downstream contraction nozzle 
was simulated by the way of a stairstep approach that represents a 45 

degree slope; it is accommodated at L/D = 2 by way of two steps, six 

cells each, in the radial direction. Even this relatively coarse 
subdivision enables useful predictions to be obtained. 

8.2 Effects of Convergence Criterion and Inlet 
Boundaries on the Prediction Accuracy 

The nature of the confined swirling flows produce the need for an 
iterative procedure for the solution of the flow governing differential 
equations, thus the accuracy of the final results depends on the con- 
vergence criterion considered. It is obvious that the smaller its value 
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the inore accurate predictions are obtained. Also, the accuracy of the 
prediction of the cited flow depends on the boundary conditions applied 
around the flow domain. Therefore, it is important to define adequately 
the boundary conditions, especially at the inlet to the flow domain. 

The effects of inlet turbulence energy (and hence the effect of dissipa- 
tion rate), inlet velocity profiles and degree of convergence are 
evaluated in this study. 

8.2.1 Effect of Refinement of Convergence Criterion 

In order to find the most economical convergence criterion that 
provides acceptably accurate predictions, a study was conducted to eval- 
uate its effect on the predicted values. Two convergence criterion were 
selected: in the former, all the normalized residual source sums for 

each variable were required to be less than 0.009; and in the latter 
were required to be less than 0.004. Results of the evaluation are 
displayed in Table XV, XVI, and XVIII. 

Parts (a) and (b) of Table XVII display the residual source sums 
for a moderate convergence criterion and a more stringent one, respec- 
tively. The flowfield axial and swirl velocity distributions, after 
achieving the displayed values, are exhibited in part (a) of Tables XV, 
XVI, -and XVIII. The former shows the flowfield predictions using part 
(a) of Table XVII while the latter was obtained from part (b) of the 
same table. Inspection of these tables reveals that the corresponding 
predicted values of the two tables are very close. The average devia- 
tions observed in the axial velocity field values for swirl vane angles 
45 and 70 degrees are approximately 1.0 percent and 1.5 percent, 
respectively. These values sharply reduce as the location moves far 
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from the centerline of the flowfield. The maximum deviation was notices 
to be about 3 percent in the centerline velocity values which are 
located very far downstream of the inlet domain. The tangential 
velocity values show a maximum deviation of less than 1 percent over the 
entire flowfield. 

Examination of the data presented in Tables XV (described in 
Section 7.5) and XVII has led to the conclusion that the predicted 
values using part (b) of Table XVII values show no preference over the 
correspondings which were obtained from part (a) of Table XVII. 

However, the former prediction method requires a computer time more than 
1.5 times that of the latter one. Therefore, to reduce the expense of 
the computer time, the convergence criterion = 0.009 was selected for 
all further predictions. 

3.2.2 Effect of Inlet Boundary Specifications 

8.2.2. 1 Effect of Inlet Kinetic Energy Profiles. To investigate 
the effect of the inlet kinetic energy k profiles on the predicted 
values, three types of specification of k were considered: 

1. Inlet k is assumed to be a constant obtained from the average 
inlet axial velocity: 

k = 0.03 u v ,o 2 (8.1) 

2. Inlet k is allowed to vary with radius and is a function of 
local inlet time-mean axial velocity at each radial location 

k = 0.03 u 0 2 (8.2) 


63 


3. Inlet k is allowed to vary with radius and is a function of 
local inlet total time-mean velocity at each radial location 

k = 0.03 (u 0 2 + v Q 2 + w 0 2 ) (8.3) 

The factor 0.03 in these equations is the recommended constant given an 
inlet relative turbulence intensity of 0.245. Predictions using these 
three inlet k assumptions with realistic inlet velocity values (measured 
u, v, and w velocities are given in part (a) of Tables VX and XVI, and 
in Tables XIX and XX. Parts (a) and (b) of Tables XIX and XX display 
the predicted axial and swirl velocities for swirl vane angles 45 and 70 
degrees, respectively. These tables illustrate that there is little 
change in velocity predictions with changes in the inlet k 
specification. The third expression for the inlet k looks like the most 
realistic one, therefore, it has been chosen for all further 
predictions. 

8. 2. 2. 2 Effect of Inlet Velocity Profiles. Numerical predictions 
of confined turbulent swirling flows are presented using various inlet 
velocity starting conditions for the case of swirl vane angles equal to 
45 and 70 degrees. In all cases considered, the inlet kinetic energy of 
turbulence k is specified as in the third case of the previous sub- 
section, that is the total local time-mean velocity magnitude has been 
used in specifying k at any radial position across the inlet domain. 

The inlet dissipation rate is specified as in a generally accepted way 
(14, 58). The validity of the flowfield predictions resulting from the 
choice of inlet profiles is assessed by comparing the predicted velocity 
profiles with corresponding experimental velocity profiles (10, 37). 
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Four possible specifications of the inlet velocities are 
considered: 

Case 1 . Flat inlet axial and swirl velocities with radial velocity 
zero are assumed, that is, both u and w are constant valued: 


u 0 f constant 
w 0 = u 0 tan <j> 

where <J> is the swirl vane angle. 

Case 2. As Case 1, except that the inlet swirl velocity profile is 
assumed to be that of solid body rotation: 


u 0 = constant 
w o = «m,o • r / R 

where w mj0 is the maximum orifice values of w which occurs at the outer 
edge r = R of the inlet. The value of w m>0 is so chosen as to maintain 
the swirl number S' the same as Case 1 for similarity purposes is now 
described. Defining 

S -iTT72 < 8 ’ 4 > 

X ' 

where 


G 0 = / (puw + pu'w' )r Z dr (8.5) 

o 

is the axial flux of swirl momentum, including the x-e direction turbu- 
lent shear stress term 


00 ■ — I , . 

G x * / (pu 2 + (p-pj) rdr 
o 


( 8 . 6 ) 


is the axial flow of axial momentum including the x direction turbulent 
normal stress term and a pressure term. 
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One useful deduction is possible when solid body rotation plug flow 
is assumed at the nozzle. That is, axial velocity u is a constant flow 
profile and swirl velocity w increases from 0 (at r = 0) to w m>0 (at r = 
d/2), the outer wall of the nozzle. Then local static pressure p and 
local swirl velocity w are related via 

P ' P» = *7 p(w2 - w mo ) 

where 

p„ = static absolute pressure at r = d/2 at the inlet station, 
x/D = 0. 

If the pressure contribution to G x is retained in the form of a 
w^/2 term, but the turbulent stress terms are omitted, the analysis 
leads immediately to 

G e ‘ 7 •’ u o “mo (d/2)3 
G x = | p u 0 2 (d/2) 2 (1 - (G/2) 2 ) 

where G * w m0 /u 0 represents the ratio of maximum velocities measured at 
the exit plane. Thus the swirl strength can be inferred from 

S = — (8.7) 

1 - (G/2 ; 

if pressure is included in the G x definition, and 

S' = G/2 (8,8) 

if pressure is omitted in G x . Both S and S' are called swirl numbers, 
but a given velocity distribution gives different S and S' values. The 
theoretical S vs. G seems to be valid only for low swirl strengths with 
S <; 0.2, see Refs. 1 and 113. For higher degrees of swirl, however, the 
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axial velocity distributions deviates from plug flow considerably - the 
major portion of the flow leaves the orifice near the outer edge (1, 3, 
11, 7,13, and 114). Then measured S valuer (including pressure contri- 
bution) are higher (114) than measured G vakniS would suggest via Eq. 
(8.7). In many practical experiments. S' values are calculated and 
quoted based on experimental swirler exit velocity profiles only, with 
the pressure contribution omitted, see for example. Ref. 115 Pressure 
variations still exist, although quoted S' values do not make use of 
these values. Under these circumstances, for a vane swirler with swirl 
vane angle <j>, it can be shown (116) approximately that 



i - (<y d > 3 
} - 


tan $ 


(8.9) 


with plug flow solid body rotation exit velocity profiles. Table XXI 
reveals how <j>. S' and G are related under such conditions for the swirl 
vane angles of special interest. 

In this study, the swirl velocity is deduced in the manner 
[$ ~ S' ~ G], with appropriate radial gradients being automatically set 
up in the Sumputer program. Notice that a gross error would be made if 
Eq. (8.7) were to be used to set up the inlet velocities [<)> ~ S ~ G], 
since Eq. (8.9) is valid for S' values, not S values. In fact, corres- 
ponding G values would be much lower than those of Table I, with G 
asympotically approaching 2 and <j> and S tend to 90 degrees and 
respectively. 

Case 3 . Measured inlet axial and swirl velocities are used (28) 
with radial velocity assumed to be zero. 
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Case 4. Measured inlet axial, radial and swirl velocity values are 


used, taking data from recent hot wire data in close vicinity of the 
swirl er exit. 

Consider first the flowfield resulting when the inlet swirl vane 
angle is 45 degrees. Figures 26a through 29a sho- predicted velocity 
profiles at various downstream axial stations, obtained when the inlet 
velocity profiles are specified by Cases 1 through 4 respectively. 
Inspection of these figures is quite revealing and may be assessed in 
the light of pitot probe and hot wire experimental data (9, 10). 

Cases 1 and 2 for swirl vane angle <|> = 45 degrees covered in 
Figures 26a and 27a give central recirculation zones terminating at 
about X/D = 2.0, with the first case having a wider and slightly longer 
zone. The central recirculation zones extend beyond 2.5 chamber 
diameters downstream of the inlet domain for the Cases 3 and 4, for the 
same swirl strength as can be seen in Figures 28a and 29a. Case 3 shows 
a thicker boundary layer at the wall than case 4. Initial spreading 
rates vary considerably: only in Cases 3 and 4 does the central recir- 

culation zone begin immediately on entry to the larger chamber, with 
Case 4 spreading most rapidly in the initial region. Cases \ and 2 do. 
not possess enough centrifugal effect because of their unrealistic inlet 
swirl velocity profiles. Also, Cases 1 through 3 do not have a radial 
component of velocity to encourage inlet radial spreading of the stream- 
lines. These inlet flow ideas may be confirmed by observing the size of 
the corner recirculation zone. Cases 2 and 4 exhibit shorter corner 
zones, with only Case 4 also possessing the correct central activity 
near the inlet -- rapid spreading with a central recirculation flow 
beginning immediately. Comparison of the predictions with the 
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experimental data (9, 10} as shown in Figure 13 clearly indicates that 
the inlet conditions of Case 4 are superior in allowing realistic flow- 
field predictions. 

Figures 26b through 29b correspond to Figures 26a through 29a but 
with the swirl vane angle increased to 70 degrees. Now the strong 
centrifugal forces present in the incoming flow play their part. 

Initial spreading rates are very high with no corner recirculation zones 
in all cases, except that of Case 3 seen in Figure 28b. All central 
recirculation zones begin immediately at the inlet and are much wider 
and longer than those with = 45 degrees. Cases 1 and 2 predicted the 
centerline velocity more closely to the experimental data than Cases 3 
and 4. However, the swirl velocity profiles are better predicted in 
Cases 3 and 4 than the other cases. 

A general deduction may be made from Figures 26 through 29 when 
compared with the experimental data (9, 10). Results demonstrate that 
the type of inlet velocity profiles shows a clearer effect on the pre- 
dictions of the 45 degree swirl vane angle flowfield than the corres- 
ponding one of 70 degrees. Also, it shows that realistic predictions 
are forthcoming only from the inclusion of the most accurate axial, 
radial and swirl velocity profiles as inlet conditions. Similar 
conclusions were reported by Abujelala and Lilley (2) using the STARPIC 
computer program along with the standard k-e turbulence model. 

8.3 General Predictions for Flowfield with D/d = 2 

The effects of swirl strength and downstream contraction nozzle are 
discussed in this section. In discussing the effects of swirl strength, 
the velocity field predictions for swirl vane angles of 0, 45, and 70 
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degrees, using measured swirler exit u, v and w velocity profiles in 
each case, are considered. Two ratios are employed in the investigation 
of the influence of downstream contraction nozzle. As mentioned 
earlier, the predicti ;s are obtained via the STARPIC computer code 
along with the modified k-e model; the model parameter formulation is 
described in Section 7.4. 

8.3.1 Effects of Swirl Strength 

The effects of swirl strength on the flow pattern are firstly 
investigated via a flow visualization technique, which uses neutral ly- 
buoyant helium-filled soap bubbles, as shown in Figure 30. This figure 
shows the swirl effects qualitatively. To obtain quantitative values, 
predictions were made for different swirl strengths. The velocity field 
predictions for swirl vane angles 0, 45, and 70 degrees are displayed in 
Figure 23 and parts (a) and (b) of Figure 29, respectively. The corres- 
ponding three-dimensional representations are given in Figure 31. Here 
the ordinate of normalized axial velocity is shown impressively as a 
function of normalized flowfield position. The flow is from right to 
left, with the sizes and shapes of corner and central recirculation 
zones being clearly evident, along with the axial and swirl velocity 
magnitudes. 

The predicted effects of swirl shown in Figures 23 and 29 confirm 
in general the well-known ideas about swirl effects on axisymmetric 
turbulent confined jet flows (115). Under nonswirling conditions a 
large corner recirculation zone exists which extends approximately to 
x/D =2.3 for the expansion geometry D/d = 2 and to x/D = 1.25 for 
D/d = 1.43, see Refs. (93) and (117). Both these results are consistent 
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with an attachment point about 8 step heights downstream of the expan- 
sion plane, as found by other researchers (9, 20, 26), The centerline 
axial velocity changes gradually from its inlet value as downstream 
development occurs. However, as the degree of inlet swirl is increased 
to 45 degrees, axial velocity profiles change dramatically. Near the 
inlet a central toroidal recirculation zone appears and the corner 
recirculation zone shortens considerably. Under strong swirl condition 
of equal to 70 degrees, a much wider central recirculation region is 
established. It promotes a very large forward velocity near the 
confining walls rather than a corner recirculation region. 

Generally, the following results can be drawn from Figures 23 and 
29: 

1. The higher the swirl strength, the higher the centrifugal 
force, hence the larger the central recirculation zone core diameter; 

2. The stronger the swirl strength, the longer the central zone 
extends downstream; 

3. The corner recirculation zone size reduces sharply as swirl 
strength increases. 

8,3.2 Effects of Downstream Contraction Nozzle 

This section discusses the prediction of the flowfields with a 
downstream contraction nozzle located at L/D = 1 and 2 for nonswirling, 
moderate and strong swirl strengths, which correspond to swirl vane 
angles = 0, 45, and 70 degrees, respectively. Two nozzles of area 
ratio 2 (weak) and 4 (strong) are being used. The weaker one has its 
upstream face contoured in a quarter circle; the stronger one has its 
upstream face in a 45 degree slope. • These are typical of exit nozzle 


71 


for ramjet and gas turbine combustors, respectively, and shown 
schematically in Figure 2. 

The nozzle area ratio and location effects on the predicted flow- 
field, using the STARPIC computer code with the standard k-e model, were 
discussed at length by Abujelala and Li 1 ley (2). They reported that the 
weak nozzle produces minor effects on the predicted flowfield, espe- 
cially at low swirl strengths and at high L/D value. The finding can be 
seen in Figure 30, which displays the flowfield photographically. 
Therefore, only the strong nozzle effects are explored numerically in 
the present study, with the nozzle located at axial stations L/D = 1 and 
2 downstream of the domain inlet.. 

Under nonswirling conditions, the blockage effects are minimal at 
either nozzle location — the central forward flow is accelerated and 
the large corner recirculation remains. Under intermediate and strong 
swirl conditions, however, the effects of the 4 to 1 area reduction are 
dramatic at either nozzle location. The corner reci rculation zone is 
encouraged to occur. This is expected due to the grid size near the 
nozzle wall; the nozzle should accelerate the flow. The central recir- 
culation zone is removed, particularly for the moderate swirl strength, 
and a small annular recirculation zone is formed for the strong swirl 
strength flowfield, as also found experimentally (9, 10), and 
numerically (48). These effects are seen in parts (a) through (c) of 
Figures 32 (with L/D = 1) and 33 (with L/D « 2) for swirl vane angles 
4> = 0, 45, and 70 degrees, respectively. These figures show axial and 
swirl velocity profiles and may be compared with the earlier Figures 23 
and 29 which were for the open-ended flowfield. The strong contractions 
also result in very large swirl velocity magnitudes and gradients in 
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central core regions, which prevail throughout the flow domain for both 
moderate and strong swirl strengths. 

Generally, the predicted nozzle effects on the flowfield character- 
istics are in good agreement with the experimental data {9, 10), except 
the values just upstream of the nozzle obstacle; the predictions show a 
small corner recirculation zona. However, the experimental data display 
an accelerated flow near the wall. It is expected that better predic- 
tions will be obtained if the grid size is further refined, especially 
near the nozzle location. 

8.4 General Predictions for Flowfield 
with D/d = 1.5 

8.4.1 Effects of Swirl Strength 

The time-mean axial and swirl velocity profiles were predicted for 
nonswirling, moderate and strong swirl flows which correspond to swirl 
vane angles of 0, 45 and 70 degrees, respectively, as shown in parts 
(a), (b) and (c) of Figure 34, Inspection of this figure reveals that 
the size the corner recirculation zone reduces and the central recir- 
culation zone builds up as swirl strength increases. These observations 
are in accord with the results in Section 8.3.1. The predicted axial 
and swirl velocity values are higher than the corresponding ones for the 
flowfield with D/d = 2 as a result of the flow area reduction. 

8.4.2 Effects of Downstream Contraction Nozzle 

In these particular predictions, a contraction nozzle of area ratio 
4 is located only at L/D = 2 downstream of the inlet domain, since 
inspection of Figures 32 and 33 shows that the nozzle produces 


73 


comparable effects at either one or two diameters downstream of the 
inlet. The predicted axial and swirl velocity profiles for swirl vane 
angles of 0, 45 and 70 degrees are displayed in parts (a), (b) and (c) 
of Figure 35, respectively. 

Comparison of the predicted velocity profiles displayed in Figure 
33 and the corresponding ones of Figure 35 reveals that the presence of 
the nozzle produces effects similar to those reported in Section 
8.3.1. However, the absolute velocity magnitudes are higher. The small 
central recirculation zone which appears in Figure 33c is not seen in 
part (c) of Figure 35, for the case of 4> = 70 degrees. 

8.5 General Predictions for Flowfield 
with D/d = 1 

8.5.1 Effects of Swirl Strength 

The previous sections describe a flowfield with an abrupt expan- 
sion. That flowfield has been found to be characterized by a complex 
flow with recirculation zones. This section discusses a flowfield that 
is found in a circular cross-sectioned duct flow without inlet expan- 
sion. The predicted results are computer plotted and displayed in 
Figure 36. Three swirl strengths vary from zero to strong swirl as 
described in Section 8.3.1 are shown in parts (a), (b) and (c) of this 
figure. Figure 36 shows no corner recirculation zone for nonswirling 
and moderately swirling flows. However, the strongly swirling flow 
($ = 70) shows a very small central reel rculation zone close to the 
domain inlet. This is expected due to the balance between the kinetic 
energy of the flow and the centrifugal force at low and moderate swirl 
strengths. The aforementioned observations, extracted from Figure 36, 
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are in accord with expectations. Vortex breakdown may occur only in 
high swirling flows (6). 

8,5.2 Effects of Downstream Contraction Nozzle 

The predictions presented in this section are similar to the ones 
of Section 8.3.1 and 8.4.1, except that there is no flow expansion in 
the present case, and more uniform flow approaches the nozzle. The 
predicted profiles for time-mean axial and swirl velocity for swirl vane 
angles of 0, 45 and 70 degrees are displayed in Figure 37. Inspection 
of this figure reveals that the general effects of the presence of the 
nozzle are again expected. The flow accelerates as the nozzle is 
approached. The central recirculation zone disappears, even at the 
highest swirl strength considered. 

8.6 General Predictions with Standard 
k-e Model 

Previous predictions shown so far have been obtained with the C p 
formulation model and general flow trends exhibited. How do these 
predictions compare with standard k-e model predictions? It is inter- 
esting to see and compare the possibilities. Indeed, for the D/d = 2 
case, predictions were given in Ref. (2) for the open-ended flow and the 
flow with a strong nozzle at L/D = 2. These may be compared with 
Figures 29 and 33, respectively , of the present study. Corresponding 
results for D/d = 1.5 are given in Figures 38 and 39, and these may be 
compared with Figures 34 and 35. similar results for D/d = 1 are given 
in Figures 40 $nd 41 and these may be compared with Figures 36 and 37. 


75 


8.7 Summary and Discussion 


8,7.1 Effects of Expansion Ratio and Swirl Strength 

Sections 8.3.1 and 8.6 have shown that the expansion ratio can 
considerably affect the aerodynamics, flame shape, heat transfer and 
stability limits, especially at high expansion ratio and high swirl 
strength. However, a few studies of the effect of levels of expansion 
upon confined swirling and nonswirling jets have been found in the 
literature. The reattachment point (the point at which the dividing 
streamline strikes to the wall) has been found of particular Interest in 
studying the effects of expansion on nonswirling flow, and thus expan- 
sion effects on the axial extent of the flow separation regions (118- 
121 ). 

In this study three expansion ratios are considered: 2, 1.5, and 

1. Computer plots of predicted streamline patterns using the C p formu- 
lation model for each expansion ratio for 0, 45 and 70 degrees of swirl 
respectively, are displayed in Figures 42 through 44. Part (a) of these 
figures represent flowfield patterns for nonswirling flow, and parts (b) 
and (c) show the streamline plot for - 45 and 70 degrees, 
respectively. 

Comparison between part (a) of Figure 42 through Figure 44 reveals 
that the length of the corner recirculation zone (CRZ) increases almost 
linearly with expansion ratio. These figures show that the CRZ extends 
to about 3, 2 and 0 times the combustion chamber diameter for D/d = 2, 
1.5 and 1, respectively. These results are in accord with the finding 
of Beck and Roschke (118), who reported that the reattachment length 
should increase linearly with both Reynolds number and the dimensionless 
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geometric ratio of step height and inlet nozzle diameter. In addition, 
the maximum reattachment length and the Reynolds number at which this 
occurs, depends on the type of inlet velocity profiles and on the flow 
geometry (119). 

Parts (b) and (c) of Figures 42 through 44 show effects of the 
expansion ratio and swirl strengths on the width and length of the 
central recirculation zone. It is clear from inspection of these 
figures that the expansion ratio alters the streamline patterns of the 
flowfield. In the case of 45 degrees swirl, it seems that as the 
expansion ratio decreases, the beginning of the central recirculation 
zone moves further downstream of the inlet to the flow domain; as the 
expansion ratio approaches one, no recirculation zone occurs in the flow 
domain. However, the effects of reducing the expansion ratio on a $ = 

70 degree flowfield are seen to be that the central recirculation zone 
sharply reduces in size, becoming very small as the expansion ratio 
reduces to the nonexpanded case. 

8.7.2 Superiority of the C y Formulation Model 

Comparison of general time-mean velocity experimental data of Refs. 
(9, 10, 122), obtained using a five-hole pitot probe and hot-wire, 
indicate that the predictions using the proposed C p formulation are 
superior to these using the standard k-e model. This is especially true 
regarding the central recirculation core diameter and length, and magni- 
tudes of axial and swirl velocities. Figures 13, 14, 23 through 25, 45 
and 46, give predictions with the C p formulation model and compare them 
with experimental data. The agreement is generally very good and better 
than would be seen were the standard k-e predictions used instead. 
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8.8 Closure 


A numerical prediction using the k-e turbulence model, with C p 
formulation obtained from recent experimental data, and different inlet 
flow boundary condition assumptions has been applied to a confined 
swirling flow. Inlet flow boundary conditions have been demonstrated to 
be extremely important in simulating a flowfield via numerical calcula- 
tions. Predictions with either flat inlet profiles, solid body rotation 
or zero radial velocity are inappropriate. 

Predictions are given for a full range of swirl strengths using 
measured inlet axial, and swirl velocity profiles. The predicted 
velocity profiles illustrate the large-scale effects of inlet swirl on 
flowfields. Predictions are included for the effects of the presence of 
a strong downstream contraction nozzle and the expansion ratio on the 
flowfield. It appears that a strong nozzle has pronounced effects on 
swirl flow cases, with discouragement of central recirculation zones, 
and forward flow in highly swirled vortex core regions further 
downstream. The expansion ratio value has large-scale effects on the 
size and location of the recirculation zones. 
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CHAPTER IX 


CONCLUSIONS AND RECOMMENDATIONS 
9.1 Conclusions 

Upon the completion of this numerical study for the two-dimensional 
axisymmetric confined turbulent swirling flow, the following conclusions 
were reached: 

1. Inlet flow boundary conditions have been demonstrated to be 
i^.tremely important in simulating a flowfield via numerical calcula- 
tions. Prediction with either flat inlet profiles, solid body rotation 
or zero radial velocity are inappropriate. 

2. The presence of a downstream contraction nozzle has shown 
considerable effects on the flowfield characteristics. It appears that 
a weak nozzle has only a minor effect on the flow. In the swirl flow 
cases, a weak nozzle leads to the discouragement of central recircu- 
lation zones with stronger vortex cores downstream possessing negative 
axial velocities. A stronger nozzle has more pronounced effects on 
swirl flow cases, and forward flow in highly swirled vortex core regions 
further downstream. 

3. The degree of swirl strength and the expansion ratio have 
strong effects on the characteristics of confined turbulent swirling 
flow. In a nonswirling flow, a large corner recirculation zone exists 
in the flowfield with an expansion ratio greater than 1. However, as 
the degree of inlet swirl increases, the size of the zone decreases and 
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a central reci rculation zone appears near the inlet. Generally, the 
size of the central zone increases with swirl strength and expansion 
ratio. 

4. Neither the standard k-e model nor its previous extensions show 
effective capability of predicting confined turbulent swirling flows. 
However, either deduced optimum values of three parameters in the model 
or the emperical C p formulation obtained via careful analysis of the 
available turbulence data can provide more acceptable accuracy in the 
prediction of these swirling flows, 

9.2 Recommendations 

Turbulence modeling has taken the attention of many investigators 
over the years. The complexity of the subject hss limited the use of 
proposed models to the flow cases that they were derived for. The range 
of the applicability of these models can be extended through simple 
modifications to their parameters. This study exemplifies the extension 
of the k-e model to confined turbulent swirling flows. 

In derivation of the proposed C y formulation, the time-mean 
velocity axial gradients were ignored because of lack of detailed data. 
As mentioned in Chapter VI, this assumption may lead to false conclu- 
sions in regions of high recirculation. Therefore, it is advised to 
evaluate the effects of inclusion of the velocity axial gradients on the 
proposed model assumptions. Also, this evaluation can be obtained 
through a careful measurement of the length scale. Finally, near wall 
turbulence data are needed to examine the model performance close to the 
wall. 
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29/1983 Axisymetric turbulent reacting 2-0 p, u, v Velocity and temperature profiles, 

and nonreacting swirling flows. k-e, Two-step and carbon monoxide and nitrogen oxide 

global reaction distributions. 
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36/1981 Coaxial jets discharging into Laser doppler Mean and fluctuating velocity and 

a large axisymetric duct. velocimeter concentration distributions. 
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TABLE III 


THE FORM OF THE COMPONENTS OF 
THE LINEARIZED SOURCE 
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A * 

S ^/Vol. 
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S */Vol. 
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S u _ i£ 
s 3x 
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->~2 
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s v + £iL _ i£ 
r 3r 

w 

y 
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ovw w 3 , s 

r ~ 2 3r 

r 

k 

li/o k 

- S.V 2 ’ 1 / 11 

G 

e 

y/o_ 

- C^pe/k 

CiC^Gpk/y 


In this table certain quantities are defined as follows: 
c u 3 , 3u. . 1 3 . 3v s 

s ' » 111 35 5 + 71? (r “ li > 


c v 3 , 3u, , 1 3 , 3v , 

s ■ aF 1 + r T? (rl1 aV > 


6 > + <!?-§> 


+ {r ^ ( ?» 2 + ( s > 2 1 


where y = ^ef c anc * stan ds for the cell control volume. 
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TABLE IV 


THE k-e EMPIRICAL CONSTANTS 
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TABLE V (Continued) 
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TABLE XI 


DISTRIBUTION OF C? VALUES OVER THE FLOWFIELD OF * = 45 DEGREES 

USING EQS. (5.4) AND (7.1) 
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TABLE XII 


DISTRIBUTION OF C 2 VALUES OVER THE FLOWFIELD OF * « 70 DEGREES 

USING EQS. (5.4) AND (7.1) 



TABLE XIII 

THE k-e OPTIMUM CONSTANTS FOR SWIRLING FLOWS 
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TABLE XIV 

PREDICTED NORMALIZED k and C VALUES USING EQ. (7.2) 
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TABLE XIV (Continued) 
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PREDICTED AXIAL AND SWIRL VELOCITIES FOR SWIRL VANE ANGLE 

* • 45 DEGREES 
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(a) Using Eq. (7.2) 


TABLE VX (Continued) 
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(b) Using Table XIII 
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PREDICTED AXIAL AND SWIRL VELOCITIES FOR SWIRL VANE ANGLE 

d. - 70 DEGREES 
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(a) Using Eq. (7.2) 


TABLE XVI (Continued) 


OjtlCOtClOClVBCOf-'C 

(Br'nonnicinf- 

nnnnn-CO-n 

dbbodooooo 


wt^noDOrrr»«P '<2 

oniio^Nttoeo 

nnmmn — 

6666666666 

* i i 


n 


nOiowMC'io^ojf* 

wvffiw^O^h^O 

wnnn* , 0 'W^ 

6666666666 

■ * i 



Oa)o<rom-- , ‘« 

Ml»n-nhOhn 

u> umi n n o o n 

6666666666 
i i i * 



tnnn 
— to tn o 

oooo 


oinontt 

t t Ifi B 31 t 
ffl co r« tfl oc O 

666666 


r>cNtnoo)<no>w«»- 

tptntn^TOcMincotin 

oooooufflipno 

d--*--ddddd 


■ 


i 


nMnOBionninffl 

aifflnininottcgin 

aoooooni>no 

00""0000 



n rt w «- cn « n co o r> O 
NtOu ji« -o-ie 3 
fflOOOOOOCDVO '*■* 

dd^^-^^odd * 


in 

O 


iS«BBaNOnc>ii» 

cof-inoO'-w^cnin 

6666666666 

i i i i i 


« 

d 


U>fflU)TTr-®in 

Of**-a®cocn®n® 

mmoiciC'-OOtO 


OOOO- — ooo 


o 

o 


8 


OOOOOPJCT®'*' 
_OOOOOCDOf«o 

ooooosatno 

OOOOOOnOOO 


o 

d 


OOOOOO'OBV 

gooooo-uio* 

OOOOOOOMf— 

ddoddddddd 



ilit- tp in mo n m m n 

01 tOl 8 «ntP-#' 

,t<tt nnwcc-oo 
dbbbbobdbb 


U“ r- W Mi r*- Mi t~J Mi Mi « 

ttvonnw-oo 

dodddddddd 


138 


(b) Using Table XIII 


TABLE XVII 


NORMALIZED RESIDUAL SOURCE SUMS 



UMON 

VMON 

WMON 

MASS 

45 

0.0090 

0.0039 

0.0086 

0.0068 

70 

0.0075 

0.0028 

0.0087 

0.0027 


(a) Moderate Convergence Criteria 


r 

UMON 

VMON 

WMON 

MASS 

45 

0.0040 

0.0016 

0.0035 

0.0034 

70 

0.0034 

0.0010 

0.0040 

0.0006 


(b) Fine Convergence Criteria 
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PREDICTED AXIAL AND SWIRL VELOCITIES USING FINE CONVERGENCE CRITERIA 
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70 deg 


PREDICTED AXIAL AND SWIRL VELOCITIES USING EQ. (8.1) 
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45 deg 


TABLE XIX (Continued) 


r*ini»TrOin<0a>iou> 

r»BT7 0!nniflf'O» 

M ra « n « — OOn« 

bbdbdbbbdb 

i i i 


«<Dd*r- — cTOf*w<i' 
n n o o rs ~ O c* a 
oddddooooo 


OBffl-OWMDOI/l 

miS'T'-wiOcooiO'T 

o)oooo)«Mnno 

d---oddodd 


o o 

— 

X — 


o 

o 


i*.O«irOinc'ii 0 O®« 

ddddddddod 

tiii 


anu)n 070 ooa<« 
OaMDi'icnQ — — mm 
ui(i»n'-0''O»i(i 

ddddddddod 

iiii 


«ncoii)inO-PiM 5 

n'Dwmi'.mwrjmm 

OMfiMO'^n^inio 

ddddddddod 

i i i i i 


> O O pi n it n 
IQQBOf'P 

» o O o» <» n o 
ddddddddod 


o 

s> 


o 

s 

a 


o 

n 


O O 

X — 


o 

d 


in n- 10 ® r* ui n m in « 
owoioMt-vf'* — 
|.»»vnn ««' , 00 

ddddddddod 


ouinn vwdmipi 
— cjcr^iMificoaltnin 
noooOBioipno 

o — — — — boob© 


>, «roou»op*t'i«® 
■ u> cm in r> n id id a in 
o ooooomi-no 

bd----bddd 


cn 

<L 

T» 


in — cmOidOcmidOP) _ 
a — a id o m « — « id ° 
mhoo-OObmo ,3 

boo- ---dob * 


«P»0®CT®<D«fi®V 
nmaaovoiiAn 

f^ajcoo — — Ocotj 

boob-- -bo 



- o omd 

— mo® 
o o n — 


o 

ar 


in*-iDLnr~inninmiM 

oi<?OiDNMrr-co<- 

VtdBOMCI-OO 

ddddddddod 


143 


0.06. 


PREDICTED AXIAL AND SWIRL VELOCITIES USING EQ. (8.2) 
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45 deg. 


TABLE XX (Continued) 
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70 deg 


TABLE XXI 


+. S' AND G CORRESPONDENCE FOR PLUG FLOW 
SOLID BODY ROTATION SWIRLERS WITH 
d n /d = 0.25 


♦ 

S' 

G 


0 

0 

0 


38 

0.547 

1.094 


45 

0.700 

1.400 


60 

1.212 

2.425 



70 


1.923 


3.846 


APPENDIX 8 


FIGURES 


o m&m: pa® ® 



Figure 1. Typical Axi symmetrical Gas Turbine Combustion Chamber 


Source t 


Anon. Tne Jet. Engine. 3rd ad. Pub. Re£. 
T.S.D. 1302 Derby, England: Rolls Royce 
Limited, 1969. 
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Figure 4. Flow Chart of the STARPIC 
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and 
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Computer Code 
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Figure 14. Predicted Velocity Profiles for Swirl Vane Angle 4 = 70 deg. 
Using Eq. (7.2) for C Variation as Compared with 
Experimental Data [ a — R e f. 9, a —Ref. 10] 
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Figure 19. Predicted and Measured Centerline Velocity Longitudinal 
Variation for Swirl Vane Angle <j> = 45 deg. 

[Predicted: Table IV; — Eq. (7.2); Measured: 

o— Ref. 9; ^ — R e f. 10T 
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Figure 21. Predicted and Measured Turbulence Parameters for Swirl 
Vane Angle <j> = 45 deg. 

IQ, *,o Table IV, — Eq. (7.2)] 




900B a 


Figure 22. Predicted and Measured Turbulence Parameters for Swirl 
Vane Angle <f> = 70 deg. 
l>, 0 --Ref. 10,*, o Table IV, Eq. (7.2)] 
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files Using Measured Inlet Axial. 
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OF POOR QUALITY 



without nozzle 



L/D = 1 



L/D = 2 

(a) $ = 0 deg. 

Figure 20. Flow Visualization of Flcwfield with D/d = 2, With and 
Without Weak Nozzle 
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without nozzle 



L/D = 1 



L/D = 2 

(b) $ = 45 deg. 

Figure 30 (Continued) 
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L/D = 2 

(c) <♦> = 70 deg. 
Figure 30 (Continued) 
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Figure 32. Predicted Velocity Profiles for Flowfie 
and Strong Nozzle at L/D = 1 







Figure 33. Predicted Velocity Profiles for Flowfield with 0/d 
and Strong Nozzle at L/D = 2 
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Figure 37 (Continued) 



ure 3b. Predicted Velocity Profiles for Flowfield with 0/d 
Using Standard k-e Model 
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and Strong Nozzle at L/D = 2 Using Standard k-e Model 













Figure 








with D/d = 1, $ = 0 deg. (No Swlrler) [ » ~R e f. 122] 
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